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I.   INTRODUCTION 

In  many  simulation  applications,  it  is  required  to 
generate  bivariate  random  variables  which  have  identical 
marginal  distribution. 

For  example,  in  reliability  problems,  the  assumption  that 
two  components  have  independent  exponential  failure  times  is 
very  often  unrealistic,  and  much  effort  has  gone  into  deriving 
bivariate  exponential  random  variables  to  handle  these  situa- 
tions [Gaver  (1972),  Olkin  &  Marshall  (1967),  Downton  (1970)]. 
Now  except  in  very  specific  physical  situations  it  may  be 
difficult  to  specify  the  complete  bivariate  distribution  of 
life  time  of  each  component.   However,  it  may  be  realistic 
to  specify  the  marginal  distributions  and  some  measure  of 
dependence   (usually  the  correlation  coefficient)  between  the 
life  time  of  each  component.   In  this  kind  of  situation,  we 
can  use  bivariate  random  vectors  having  given  marginal  dis- 
tribution and  dependence  to  solve  the  problem  in  simulation. 
To  generate  these  vectors,  there  exists  some  previous  works 
and  existing  methods,  discussed  in  Section  II.   As  seen  in 
Section  II,  most  previous  work  is  specific  to  specified 
marginal  distributions  and  uses  inverse  transformation  methods 
as  a  basic  concept. 

An  example  is  the  recent  work  by  Johsnon  and  Tenenbein 
(1979) ,  to  generate  a  bivariate  random  vector  (X,Y)  which  has 
marginal  distribution  F,(x),  F2(y)  and  correlation  p ,  by  a 
weighted  linear  combination  method. 


Define 


X   =   F^-^(H^(U)) 


Y   =   F2-^(H2(V)) 


or 


Y   =   F^-'-d  -  H^CV)) 


where  H,  and  H^  are  the  cumulative  distribution  functions 


(c.d.f.)  of  U  and  V  respectively  and 


U   =   U' 


V   =   cU'  +  (1  -  c)V' 


where  U',  V  are  i.i.d.  random  variables  with  probability 
density  functino  g(  ).   In  this  procedure,  F,  ,  F„  ,  and  c 
are  specific  to  the  marginal  distribution  and  correlation 
desired.   The  functions  F.   and  F^   are  difficult  to  compute 
in  most  cases  and  the  weighting  factor  c  is  also  difficult  to 
calculate.   Moreover  most  of  the  work  in  univariate  random 
number  generation  has  been  aimed  at  avoiding  having  to  calcu- 
late inverse  cumulative  distribution  functions  such  as 
F,  (  .)  and  F-  (•).   These  are  the  reasons  why  many  proposed 
methods  are  specific  to  a  specified  marginal  distribution. 


Again  special  properties  of  certain  random  variables  such 
as  infinite  divisibility  have  been  exploited  to  give  easily 
generated  bivariate  random  variables,  often  though  with 
limited  ranges  of  dependency.   One  very  clever  scheme  by 
Gaver  (19  72)  to  generate  bivariate  exponential  random  varia- 
bles uses  the  fact  that  the  sum  of  a  geometrically  distributed 
number  of  exponential  random  variables  (Y)  is  exponentially 
distributed  and  that  the  minimum  of  this  geometrically  dis- 
tributed number  of  independent  logistic  random  variables  (Z) 
is  exponential.   Clearly  when  Y  is  large,  Z  is  small.   This 
scheme  is  of  course  very  specific  to  exponential  marginal 
distributions  and,  via  an  exponential  transformation,  to 
uniform  random  variables.   To  avoid  these  kinds  of  limitations 
and  to  make  the  generation  of  bivariate  random  variables 
simpler  and  more  automatic  in  simulations  we  develop  here  a 
scheme  presented  in  Jacobs  and  Lewis  (1977) . 

This  scheme,  the  mixture-truncation  method,  is  a  very 
general  tool  which  requires  only  that  a  method  be  available 
for  generating  random  variables  with  the  desired  marginal 
distribution.   The  mixture-truncation  method  scheme  for 
generating  bivariate  random  variables  is  as  follows. 

Let  F(x)  be  the  common  marginal  distribution,  of  the 
bivariate  random  variable  (Y,Z) .   Let  p  be  the  desired 
correlation  between  Y  and  Z  (which  may  or  may  not  be  attain- 
able generally  or  in  particular  with  the  mixture-truncation 
scheme) . 


Define  the  transition  matrix  P  as 


P   = 


1  -  a. 


l-a2     o.^ 


with  stationary  vector 


IT   P    =    TT 


1  -  a, 


1  -  a. 


1-a,  +l-ap'  1-a,  +l-a2' 


and  let  the  range  of  the  random  variable  X  with  distribution 
function  F(x)  be  x*   Then  generate  (Y/Z)  as  follows. 

1.  ( Initiallization) 

i)   Choose  an  "allowable"  x   from  (x,,x  )  (the 

o  '<,   u 

allowable  range  of  x  will  usually  be  smaller 
than  X  and  depends  on  p  and  F(x)) . 
ii)   Set  7T,  =  F(x  )  ,  tt_  =  1  -  tt,  and  compute  a,  and 

iii)   Denote  by  X,  the  random  variable  X  truncated  to  the 
left  of  X   and  by  X  truncated  to  the  right  of 

X  . 

o 

2.  (Generate  Y)  .   Choose  Y  from  X,  with  probability  tt, 
or  choose  Y  from  X^   with  probability  7t_. 

3.  (Generate  Z) 

i)   If  Y  is  chosen  from  X,,  choose  Z  from  X,  with 

probability  a   or  choose  Z  from  ^^   with  probability 
1  -  a 
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ii)   If  Y  is  chosen  from  X^ ,  choose  Z  from  X^  with 
probability  1  -  ct2  or  choose  Z  from  X^  with 
probability  o.^- 

In  Section  III,  we  will  show  that  the  correlation  p  between 

Y  and  Z,  if  X   is  properly  chosen,  is 


P 


=   E[[^-^^^>]  [ZJLMZI] 


a^         02 


=   6  M, 


where 


=   a-,  -  (1  -  a-)  / 


M   = 


JM^   -  1^2^   ^1  ^2 
2 


n^   =   E{X^},    ^2      =      E{X2}  , 


a^^   =   VAR[X]  , 

Moreover,  the  generated  bivariate  random  vector  (Y,Z)  will 
have  marginal  distributions  F(x)  and  correlation  p  and  the 
joint  distribution  function  of  (Y, Z)  will  be 
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F(y,z)      = 


F(z)    F(y) 


if      z  <  X    ,    y  <  X 
—    o      -^  —    o 


(1-a    ) 
{a,    +   ^[F(z)-F(x^)  ]  }F(y)       if      z  >  x^,    Y  <  ^^ 


(1-a    ) 
{a,    +   ^[F(y)-F(x^)  ]  }F(z)     if      z<x    ,    y>x 

±  TT  ^  O  —      O  O 


a^TT,  +  (1-a    )  [F(z)-F(x    )  ]    --        if      z<x    ,    y>x 

11  1  O  7T«  O  O 


+    {(1-a^)  +— [F(z)-F(x^)]}[F(y)-F(x^)] 

2  TT-^  O  O 


These  relationships  will  be  developed  in  Section  III.   Note 
that  (Y,Z)  may  be  continuous  or  discrete  random  variables 
or  mixtures  of  both,  though  in  this  thesis  we  concentrate  on 
continuous  cases. 

The  key  problem  in  this  very  simple  algorithm  comes  at 
the  initialization  steps  i)  and  ii) .   There  are  two  degrees 
of  freedom  in  the  selection  matrix  P  but  setting  tt,  =  F(x  ) 
constraints  reduce  to  one  degree  of  freedom.   Specifying 
a  desired  correlation  further  constrains  the  degree  of  free- 
dom, though  not  completely.   Subject  to  the  constraint  that 
a,  and  a^  are  probabilities  there  may  be 

a.  no  values  x  which  will  give  (Y,Z)  with  correlation 

P 

b.  one  value  x 

o 

c.  a  range  x  • 
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The  main  numerical  problem  of  initialization  step  (i)  then 

is  to  compute  the  range  of  allowable  x   for  a  given  p  and 

F(x)  . 

The  main  statistical  problem  is  then  to  choose  which  of 

the  bivariate  random  vectors  (Y,Z)  indexed  by  x   e  x  to  use. 

An  alternate  solution  to  the  statistical  problem,  giving 

another  algorithm  is  then  to  let  x  have  some  distribution  in 

o 

the  range  x^r   possibly  uniform  or  triangular,   this  not  only 

alleviates  the  problem  of  picking  a  particular  x   e  x   ^^"^ 

it  also  smoothes  out  the  distribution  of  (Y,Z)  and  possibly 

makes  it  continuous.   The  computation  of  Y  fo^  given  o  is 

o 

illustrated  in  subsequent  sections  for  uniform,  exponential 
and  gamma  marginals  for  (Y,Z). 

Before  doing  this  and  developing,  in  Section  III,  the 
results  already  given   here,  we  review  a  few  existing 
methods  for  generating  bivariate  random  variables  in  Section 
II.   These  are  methods  which  seem  fairly  tractable,  for  use 
on  a  computer. 
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II.   REVIEW  OF  EXISTING  METHODS 

There  are  several  existing  methods  for  generating 
bivariate  random  vectors,  but  most  of  them  are  specific  to 
particular  marginal  distribution  and  use  inverse  transfor- 
mation method.   The  inverse  transformation  method  is  a  very 
useful  univariate  procedure  which,  unfortunately,  is  not 
possible  to  use  with  many  distributions  because  it  is  diffi- 
cult and/or  uneconomic  to  compute  the  inverse  functions.   We 
survey  here  some  of  the  methods  which  are  germane  to  this 
thesis,  in  particular  concentrating  on  bivariate  random 
variables  with  uniform,  exponential  and  gamma  marginals.   First 
we  review  the  problem  of  determining  the  range  of  correlation 
coefficient  p  which  can  be  obtained  for  bivariate  distribution 
with  given  marginal  distributions,  again  considering  only 
the  continuous  case. 

A.   RANGE  OF  CORRELATION  COEFFICIENT  p 

Suppose  that  Y,  Z  are  random  variables  with  an  arbitrary 
joint  distribution  F(y,2)  with  finite  second  moments.   Then 
in  general  the  correlation  coefficient  p  can  take  all  values 
in  the  closed  interval  [-1,1].   But  with  specific  marginal 
distribution  Fj^(y)  and  F^Cz),  the  class  of  all  F(y,z)  need 
not  attain   the   values  of  -1,  1,  of  p.   The  necessary  and 
sufficient  conditions  that  there  exist  determination  of  F(y,z) 
with  p  equal  to  1  and  -1  are  given  by  Moran  (1967)  as  follows. 
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i)   there  exist  constants  a  and  6  such  that  aY  +  6  has 
the  same  distribution  as  Z, 
ii)   the  distribution  of  Z  is  symmetrical  about  its 
mean. 
To  see  this,  rescale  Y  to  have  the  same  mean  and  variance 

as  Z .   If  there  exists  an  F(y,z)  such  that  p  =  1  we  have 

2  .  . 

E[Z-Y]   =0,  so  that  Y  =  Z  with  probability  one.   Then  Y 

must  have  the  same  distribution  as  Z.   On  the  other  hand 

if  there  exists  an  F(y,z)  such  that  p  =  -1  we  shall  have 

2 
E[Z  +  Y]   =  0,  and  Y  =  -Z  with  probability  one.   Thus  if  both 

bounds  are  attainable  Z  must  have  a  symmetric  distribution. 

Given  general  marginal  distributions  F, (y)  and  F^(z)/ 

Mardia  (1970)  showed  Frechet  bounds  as 


max[0,Fj_(y)+F2(z)-l]  <_F(y,z)  £  min  [F^  (y)  ,  F^  (z)  ]     (II-A-1) 

From  this  we  can  find  the  range  of  possible  values  of  p . 
For  simplicity  we  now  confine  our  consideration  to  distribu- 
tions F(y,z)  of  positive  random  variables  whose  deriva- 
tives F '  (y)  and  F'(z)  are  strictly  positive  for  y  >  0,  z  >  0, 
respectively.   Suppose  also  that  the  variates  are  scaled  to 
have  unit  variances.   Let  G, (u) ,  G^ (v)  be  the  inverse  func- 
tions of  Fj^(y),  F^Cz),  i.e. 


F^[G^(w)]   =   w 


where    0<w<l,    i=l,    2.      Then   the   correlation   coefficient 
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between  Y  and  Z  is  given  by  the  equation 


00  oo 


p 


I         /    y    z    dF(y,z)     -    E[Y]    E[Z]  (II-A-2a) 

0         0 


1         1 

./         /    G, (u)    G    (v)     dK(u,v)     -    E[Y]    E[Z]  (II-A-2b: 

0         0-^ 


where  K(u,v)  is  the  joint  distribution  of  the  quantities 
U  =  F, (y) ,  V  =  F^(z).   Then  U,V  are  jointly  distributed  on 
the  square  Oj<U_<_l,  0<_V_<1,  in  such  a  way  that  the 
marginal  distributions  are  unifon:n  on  the  unit  intervals. 
From  expression  (II-A-1)  the  minimum  correlation  is  attained 
when  the  probability  is  concentrated  uniformly  on  the  line 
U  +  V  =  1.   The  minimum  value  of  p  then  is 


1 

Pm-ir,   =    I  ^1  (^)  G^(l-u)  du  -  E[Z]  E[Y]       (II-A-3) 
mm      ^■'       1  2 


The  corresponding  F(y,z)  will  be  a  singular  distribution  with 
all  the  probability  concentrated  on  the  line  F, (y)  +  F^(z)  =  1 
In  fact  Y  and  Z  are  an  antithetic  pair.   The  maximum  value 
of  p  is  attained  when  the  probability  is  concentrated  uni- 
formly on  the  line  U  =  V.   Then  the  maximum  value  of  p  is 


1 


p 
max 


I   Gj^(u)    G2(u)    du   -    E[Z]    E[Y],  (II-A-4) 
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with  the  corresponding  probability  concentrated  on  the  line 

By  using  this  result  we  will  get  lower  and  upper  bounds 

of  the  correlation  coefficient  for  uniform,  exponential  and 

gamma  marginal  cases.   In  the  uniform  marginal  distribution 

case,  we  can  see  p  .   =  -1  and  p     =1  from  Moran ' s  condition 

mm  max 

ii) ,  i.e./  uniform  distribution  is  symmetric  about  its  mean. 
For  the  exponential  marginal  distribution  case  suppose  Y  and 
Z  have  exponential  distributions  with  unit  mean  and  variance. 
Then 


F^(y)   =   1  -  e"^,   F^Cz)   =   1  -  e"^ 


and  the  inverse  functions  are 


G^(u)   =   -  ln(l  -  u) , 


G2(v)   =   -  ln(l  -  V)  . 


From  equations  (II-A-3)  and  (II-A-4) 


1 
^min   "    ^    ^1^^^  G2(l-x)  dx  -  E[Y]  E[Z] 


1  2 

f  In  X  In  (1-x)  dx  -  1   =   1  -  V 
0  ^ 


-0.64493 
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max 


I  G^(x)  G^ix)    dx  -  E[Y]  E[Z] 


/  In  X  In  X  dx  -  1   =   1 


The  p  .   can  be  attained  when  all  Y  and  Z  are  concentrated 
mm 

on  the  line  e  -^  +  e   =1.   Also  the  p     can  be  attained 

max 

when  all  Y  and  Z  are  concentrated  on  the  line  e  -''^  =  e  , 
i.e. ,  Y  =  Z . 

For  the  gamma  marginal  distribution  case,  suppose  that 
Y  and  Z  have  gamma  type  distributions  with  probability  density 
functions 


fi(y)  = 


r(a) 


-y   a-1 


a  >  0   ,  y  >  0 


f2(z) 


1   -z  e-1 

e   z 


r(B) 


>  0   ,  z  >  0 


Then 


E[Y]   =   a,    E[Z]   =   3/   VAR[Y]   =   a   VAR[Z] 


The  P  •   will  be  attained  when  the  probability  density  is 
concentrated  on  the  line 


1    c^_-x   a-1 


r(a) 


ex    dx  + 


r(6)  n^ 


'  -X   3-1  -,       T 
e    X     dx   =   1 
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This  defines  y  uniquely  as  a  function  of  z  which  can  be 

written  y  =  A(z) .   The  p     is  then,  on  rescaling  the 
■^  mm 

covariance 


_   _  .„.  „ 1/2 

mm      Q 


=   {  J  u  A(u)  fj_(u)  du  -  a3}/(a6) 


When  a, 6   become  large,  p  .   tends  to  -1.   And  the  p 

mm  max 

will  be  attained  when  the  probability  density  is  concentrated 
on  the  line 


1  C         ~X         Ct~l  1  f         ~X         C6~l 

ruT  o''    ^        ^  ax     =     p-jgj-  ^/  e        X  dx 


This  defines  y  as  a  function  of  z  which  can  be  written 

y  =  B(z).   The  p     then  is,  on  rescaling  the  covariance 

■^  max  ^ 


Pt.=  ^   =   W  ^  B(u)  f,(u)  du  -  aB}/(a6)-^^^ 
max       «. '  J. 


when  a  =  3,  p     tends  to  1.   Schmeiser  and  Ram  Lai  (19  79) 
max 

showed  the  obtainable  correlations  between  random  variables 

Y  and  Z  having  gamma  marginal  distribution  with  density  function 

a.-l 
fi(x)   =   [(x/e^)  "•    exp(-x/3i)/(6ir(a^)] 


for  X  >  0,  a^  >  0,  6  .  >  0,   i  =  1,  2. 
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The  Figures  (Il-a)  and  (Il-b)  show  the  obtainable 
correlation  as  a  function  of  a.2t    given  a-j^  =  1  and  5,  for 


,  .  =  1. 

1 


q:i=5 


rrfimfiiTlil'milWMiri   Tii1»  WVlil  l^«ll'fa  irVlift-i'flJt  lIT    i    ftftMi  .1  .a.lsl 


30 


00 


ai 


Figure  Il-a.   Obtainable  correlation  as  a  function  of 


a   for  a,  =  1,  3-,  =  1/  by  Schmeiser 
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q:i=i 


«■  '•'•V  "I 


L... 


30 


CO 
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Figure  Il-b.   Obtainable  correlation  as  a  function  of  a, 


for  a,  =  5,  6-  =  1,  by  Schmeiser 
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B.   GENERAL  METHOD  REVIEW 

1.   Johnson  and  Tenenbein's  General  Method 

A  general  method  of  constructing  a  bivariate  dis- 
tribution, whose  marginal  distribution  functions  are  F, (x) 
and  F_(y),  is  proposed  by  Nataf  (1962),  can  be  represented 
as  follows. 

General  Method 

i)   Consider  any  two  continuous  random  variables  U 
and  V  with  probability  density  functin  h(u,v). 
ii)   Let  X'  =  H^(u)  and  Y'  =  H2 (v) ,  where  H^(u)  and 
H2(v)  are  the  cumulative  distribution  functions 
of  U  and  V,  respectively, 
iii)   Define 


X   =   F^^(X')   =   F^^[Hj^(u)]        (II-B-1) 


and 


Y   =   F^-'-CY')   =   F2^[H2(v)]        (II-B-2a) 


or 


Y   =   F^'^d-Y')   =   F^-^Ll  -  H^Cv)  ]  (II-B-2b) 
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since  X',  Y'  and  1-Y'  are  uniformly  distributed  over  the 
interval  [0,1],  X  defined  by  expression  (II-B-1)  and  Y 
defined  either  by  (II-B-2a)  or  (II-B-2b)  will  have  a  joint 
distribution  whose  marginal  distribution  functions  are  F-,  (x) 
and  F^ (y) .   The  method  is  probably  what  one  would  think  of 
first  but  again  involves  inverses.   Based  on  this  general 
method/  Johnson  and  Tenenbein  (19  79)  developed  two  proce- 
dures for  generating  (and  also  constructing)  bivariate  dis- 
tributions whose  marginal  distributions  and  measures  of 
dependence,  as  given  by  Kendall's  T  or  the  grade  correlation 
coefficient  p  ,  are  specified.   Note  that  these  Kendall's  T 
and  grade  correlation  coefficient  p   are  not  the  same  as  the 
ordinary  product  moment  correlation  coefficient  p  we  have 
been  considering;  these  measures  are  discussed  by  Kendall 
(1962)  and  Kruskal  (1958)  and  can  be  defined  as  follows.   Let 
X  and  Y  be  continuous  random  variables  having  some  joint 
probability  density  function.   Let  (X,  ,Y,),  (X^/Y^)  and 
(X^,Y^)  be  three  independent  pairs  of  observations  having 
the  same  joint  density  function,  then 


7   =   2  P[(X^-X2)  (Y3_-Y2)  >  0]  -  1 


Pg   =   6  P[(Xj_-X2)  (Y^-Y^)  >  0]  -  3 

The  first  procedure  for  generating  bivariate  pairs,  called 
the  ViLC    (weighted  linear  combination)  ,  defines 
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U   =   U'  (II-B-3) 


and 


V   =   cU'  +  (l-c)V'  (II-B-4) 

where  0  <_  c  <_  1,  U'  and  V  are  independent  and  identically 
distributed  random  variables  with  probability  density  func- 
tion g(')  .   Then  X  and  Y  are  obtained  from  the  general 
method  using  equation  (II-B-1)  and  (II-B-2a)  if  the  depen- 
dence measure  is  positive  or  equations  (II-B-1)  and  (II-B-2b) 
if  the  dependence  measure  is  negative. 

In  order  to  apply  the  WLC  procedure,  we  must  obtain 
expressions  for  H, (u) ,  U^iv) ,    o    (u,v)  and  r(u,v),  in  terms 
of  c  and  g(')  .   The  values  of  H,  (u)  and  Hp(v)  allow  us  to 
apply  the  general  method  for  a  given  choice  of  c  and  g(t) . 
The  expressions  for  T(u,v)  and  a    (U/V)  allow  us  to  specify 
c  for  a  given  choice  of  g ( • )  in  terms  of  the  required  value 
of  either  T  or  a  .   From  equations  (II-B-3)  and  (II-B-4) 

u 
H^(u)   =     /   g(t)  dt  (II-B-5) 


H  (v)  =      j    !   g(u')  g(v')  du'dv'    (II-B-6) 
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where 


R^   =   {(u',v'):   cu'+(l-c)v'  ^  v} 


and  the  joint  density  function  of  u  and  v  is 


h{u,v)       =      j^   g(u)  g(^x^^  (II-B-7) 


00  00 


p  (u,v)   =   12    I     I  H^(u)H2 (v)h(u,v)dudv  -  3    (II-B-8) 

—  OO        —00 


0    , 
T{u,v)       =   4    I  G„(^^-^t)  g  (t)  dt  r  (II-B-9) 


2^  c 


where 


g^Ct)   =     /  g(w  +  t)  g(w)  dw 


G^(t)   =     /  g, (X)  dx 


2 


Using  equations  (II-B-5) ,  (II-B-6) ,  (II-B-7),  (II-B-8)  and 
(II-B-9),  we  have  to  evaluate  H-,(u),  H^  (v)  ,  h(u,v),  T(u,v), 
and  p  (u.v)  for  all  values  of  c  and  for  specified  g(*). 
As  Johnson  and  Tenenbein  mentioned  in  their  report,  these 
integrals  are  quite  tedious  to  perform,  so  they  gave  some 
computation  results  in  their  report. 
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The  second  procedure,  called  TVR  (trivariate  reduction) 
is  discussed  by  Mardia  (1970) .   In  this  case  U  and  V  are 
defined  as 

U   =   U'  +  ez'  (II-B-10) 

V   =   V  +  BZ'  (II-B-11) 

for  0  <_  3  <  °°f  and  U',  V  ,  and  Z'  are  independent  and  iden- 
tically distributed  random  variables  with  probability  density 
function  g ( • ) . 

Like  the  WLC  procedure,  one  needs  to  define  H, (u) , 
H^ (v) ,  h(u,v),  p  (u,v)  and  T(u,v)  as  a  function  of  6  and  g('). 
From  equations  (II-B-10)  and  (II-B-11)  it  follows  that 


H^  (u)   =   H  (v)   =   /  /g(u')  g(z')  du '  dz '    (II-B-12) 

R. 


where 


R^   =   {  (u'  ,z  ')  :   u'  +  Bz'  <_  u}  , 


h(u,v)   =     /  g(u  -  Bz)g(v  -  Bz)  g(z)  dz       (II-B-13) 


and 
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T(u,v)   =   4    /  [G2(3t)]^  q^it)     dt  -  1       (II-B-14) 


where 


92  (t)   =     /  g(w+  t)  g(w)  dw 


G  (t)   =     /  g  (x)  dx 


2 


using  equations  (II-B-12) ,  (II-B-13) ,  (II-B-14)  and  (II-B-8) 
the  same  computations  are  needed. 

In  this  general  WLC  and  TVR  methods,  there  are  two 
problems.   One  is  the  need  for  inverse  transformations.   The 
other  is  the  need  to  generate  coupled  random  variables  to 
create  dependence  after  the  inverse  transformation  is  applied 
Doing  this  by  linear  combinations  is  not  necessarily  the 
simplest  and  most  felicitous  with  respect  to  calculation  of 
correlations  and  range  of  correlations. 

A  significant  simplification  can  be  achieved  by 
obtaining  a  "smooth"  bivariate  uniform  pair  (U,V)  whose 
correlation  can  be  varied  between  the  limit  of  correlation 
which  can  be  obtained  for  uniform  marginals.   These  limits 
are  (-1,1)  since  the  antithetic  pair  (U,l-U)  has  correlation 
-1. 
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Note  that  the  problem  is  essentially  one  of  obtaining 
a  bivariate  uniform  random  variable  with  positive  correlation 
because 


corr(U,V)   =   p 

'        ^u,v 


■Pu,l-V  =   corr(U,l-V) 


L  -p^_u  ^  =   corr(l-U,V) 


Two  fairly  simple  schemes  for  obtaining  bivariate  uniform 
distributions  are  discussed  in  the  next  section. 

C.   REVIEW  OF  SOME  SPECIFIC  METHODS 

The  mixture-truncation  method  is  a  general  algorithm 
for  generating  bivariate  pairs  (Y,Z)  with  given  marginal 
distributions  and  requires  only  the  availability  of  a  generator 
for  the  marginal  distribution  and  calculation  of  constants. 
Unless  randomization  is  used,  however,  it  does  give  a  bi- 
variate distribution  with  a  discontinuous  joint  distribution, 
as  in  Section  (III-D) .   This  could  be  a  disadvantage.   It 
is  interesting  therefore  to  compare  it  to  other  generations 
which  are  specific  to  the  given  marginal  distribution  and 
we  do  this  here  for  the  cases  for  which  the  mixture-truncation 
method  is  illustrated  in  this  thesis,  namely  uniform,  exponen- 
tial and  gamma. 

There  are  many  schemes  available  for  uniform  and  exponen- 
tial cases  though  it  should  be  noted  that  smooth,  simply 
generated  pairs  with  easily  calculated  correlations  have 
only  recently  been  available.   In  particular  we  give  here  two 
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very  recent  schemes  for  smooth  bivariate  uniforms  which  are 
easily  generated  and  whose  correlations  can  be  calculated; 
similarly  two  schemes  for  exponential  which  are  essentially 
the  uniform  schemes  after  transformation.   When  one  comes 
to  gamma  marginals  one  is  again  in  difficulty  and  we  describe 
a  recent  scheme  by  Schmeiser  and  Ram  Lai  which  involves 
however,  computation  of  the  inverse  gamma  distribution 
function  and  very  difficult  initialization  computations. 
The  mixture-truncation  method  is  quite  competitive  here.   In 
fact  if  one  can  compute  the  inverse  gamma  distribution  then 
bivariate  gamma's  can  be  computed  as  [Y  =  F   (U) /  Z  =  F   (V)] 
where  U  and  V  are  a  bivariate  uniform  generated  by  the  schemes 
mentioned  above;  this  would  be  simpler  for  any  marginal  dis- 
tributions than  some  of  the  general  schemes  mentioned  in  the 
previous  sections.  Comparisons  of  the  generating  schemes  given 
in  this  section  with  the  mixture-truncation  method  are  given 
in  Section  IV-D  and  Section  V-D. 

1.   Lawrance  and  Lewis's  Bivariate  Uniform 

Lawrance  and  Lewis  (1979)  show  that  a  simple  trans- 
formation of  the  NEAR(l)  (New  Exponential  Autoregressive) 
process  gives  a  two-parameter  family  of  Markovian  random 
variables  with  uniform  marginal  distributions .   This  method 
generates  a  correlated  uniform  pair  as  a  multiplicative  mix- 
ture of  uniform  variables,  where  the  parameters  a  and  B  take 
on  values  between  zero  and  one,  with  the  condition  that  they 
not  both  be  one.   Let  Y  be  a  uniform  [0,1]  random  variable, 
and  define 


29 


where 


Z   = 


r  e  Y' 


wp   a 


wp   1-a 


r  U 


•-  U 


(1-a) 6 


wp 


wp 


1-e 

1  -  (1-a) 


ag 
1-  (1-a) 


where  U  is  an  independent  identically  distributed  uniform 
(0,1)  random  variable.   The  correlation  structure  is  defined 
as  follows  as  a  function  of  a  and  8; 


aB 


Y,Z 


2  +  B   1  +  (1-a) 


The  model  can  be  reduced  to  a  one-parameter  model  by  suitably 
relating  a  and  Q,    e.g.,  a.   =  3  and  all  positive  values  of  p 
can  be  obtained.   Consequently  Y  and  1-Z  will  have  a  full 
range  of  negative  correlations.   A  generating  procedure  for 
this  bivariate  pair  of  uniform  random  variables  is  as  follows. 


1. 

2. 
3. 


Generating  Procedure 
(Initialization) .   For  given  correlation,  find  suitable 

parameter  values  a  and  3. 

1-3 


Generate  U,  ,  set  Y  -^  U,  ;  r  =  1-a,  P  -^  %  ,n       . 
J.  1  ±—  ^1-  a; 

Generate  \J^. 
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^2 

4.  If  U2  ^  r,  set  U  -*-  —  and  go  to  8. 

5.  Otherwise  U  ^  -^ . 

1  -  r 

6.  If  U  _^  P,  set  Z  ^  Y^  Y^  and  go  to  10. 

7.  Otherwise,  set  Z  ^  Y^(j-^)^^  and  go  to  10. 

8.  If  U  <_  P,  set  Z  ^  -  and  go  to  10. 

9.  Otherwise,  set  Z  ^  (?  I  p)^^- 

10.   Go  to  2  and  repeat  until  the  desired  number  of  bi- 

variate  uniform  variables  (Y,Z)  are  obtained. 

The  program  is  listed  in  the  appendix  and  the  scatter  plot 

and  bivariate  histogram  of  resulting  output  are  in  Section 

IV-D. 

2 .   Bivariate  Uniform  By  Transformation  of  Gaver ' s 
Bivariate  Exponential 

In  Gaver ' s  bivariate  exponential  to  be  discussed  in 

Section  II-C-4,  we  can  define  a  bivariate  exponential  random 

vector,  (Y,Z)  with  correlation  p  =  -  P/2  as  follows: 


^  =  1^^-17^ I?p^        (ii-c-i) 


and  Z  is,  conditional  upon  N  =  n,  a  gamma  random  variable 
with  shape  parameter  n  and  mean  n(l-p),  where  N  is  a  geometric 

random  variables  with  probability  density  function 

x-1 
f(x)  =  p(l-p)    ,  X  =  1,  2,  ...  and  U  is  a  uniform  random 

variable  over  [0,1].   From  this  we  can  generate  a  pair  of 
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uniform  random  variables  (V,W)  by  transformation.   Since  we 
know  that  to  generate  an  exponential  random  variable  X  from 
a  uniform  random  variable  U[0,1),  we  use  the  following 
transformation: 

X  =   -  ln(l  -  U)  (II-C-2) 

The  resulting  X  has  the  exponential  distribution  with  unit 
mean.   From  equations  (II-C-1)  and  (II-C-2) ,  we  can  generate 
V,  a  uniform  random  variable  over  [0,1]: 


^  .   (I-P)  u^/^ 

1  -  P  ui/» 


Further  also  we  can  generate  W,  a  uniform  random  variable 
over  [0,1] : 


N    1-P 
w  =  (  n  u.) 

i=i  ^ 


This  follows  since  we  know  that  Z  has,  conditionally  upon 
N  =  n,  the  gamma  distribution  with  shape  parameter  n  and 
mean  is  n(l-p) . 

In  general  the  resulting  V  and  W  will  be  negatively 
correlated  and  (V,l-W)  will  be  a  positively  correlated  pair, 
Because  the  correlation  structure  will  be  changed  by  trans- 
formation, the  resulting  (V,W)  need  not  have  the  same 
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correlation  with  (y,Z).   The  correlation  coefficient  of 
(V,W)  will  be  computed  as  a  function  of  P  as  follows. 

p   =   12  E[VW]  -  3 

where 


E[VW]   =   E  [E[VW|N  =  n]] 


1  ,,_^,,,l/n         1 


r  f       (1-p) U  ^   ,     r    (      1-P  ,  ,n, 
[  J   -^ ^  -^/^   du  -  [  J  u    du]  ] 


^^'o'   1-P  U^/^'       0 


Unfortunately  this  computation  of  p  as  a  function  of  P  is 
difficult.   As  an  example  we  used  here  the  same  function  for 
p  as  holds  in  the  exponential  case. 

Generating  Procedure 

1.  (Initialization) 

i)   Compute  P  for  given  correlation  p 
ii)   Choose  N  from  a  geometric  distribution  with 
parameter  1-P 

2.  Generate  a  uniform  [0,1]  random  variable  U,  and  define 

(1-P)  u/^" 

V  =  Vt- 

1  -  p  U3_i/^ 

3.  Generate  N  =  n  uniform  [0,1]  random  variables  U., 
i  =  1,  ...,  n,  and  define 
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n     , 

w  =  (  n  u.)   P 

i=l  ^ 


4.   Deliver  (V,W)  and  go  to  2  until  a  sufficient  number 

of  bivariate  pairs  have  been  generated. 

The  program  is  listed  in  the  appendix  and  the  scatter  plot 

and  bivariate  histogram  of  resulting  random  vectors  are  in 

Section  IV-D. 

3.   Marshall  and  Olkin's  Positively  Correlated 
Bivariate  Exponential 

Suppose  X  is  the  age,  or  length  of  service  of  the 
first  device  at  the  time  of  death,  governed  by  two  indepen- 
dent Poisson  processes  Z, (t) ,  Z,^(t)  with  parameters  X,, 
A.^  respectively  and  Y  the  same  of  the  second  device, 
governed  by  two  independent  Poisson  processes  Z_(t),  Z^-(t) 
with  parameters  X_,  A,^/  respectively.   Further  suppose  the 
second  device  is  placed  in  operating  after  a  time  6  later 
than  the  first  device.   Then  the  joint  distribution  of  (X,Y) 
is  defined  as 

P[X>x,Y>y]  =      F(x,y) 

exp{-A,  x-A2y  -  X,  ^rnax  [x,y+min  (x,  6)  ]  }         if      6    >_  0 

L  exp{-A,x-A2y-A,  2^^^  ty '^■'"^i"  (Y' <5)  ]  }  if       6     <_   0 

(II-C-3) 
If  6  =  0,  then  equation   (II-C-3)  reduces  to 
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P[X>x,Y>y]       E      F(x,y) 

=      exp{-A,x  -    A_y  -    X,2   inax(X/y)} 
Note   that   since   the  marginal   distributions    do  not   depend  on 


E[X]       =  ^  , 


VAR[X]       =      '^ J    ' 


(A3_+A^2) 


E[Y]       =  ^ 


^2    "^^12' 


VAR[Y]       = 


(A2+Xi2^^ 


We   also  have,    for    5    >_  0 

-(A^+A,2) 

1  ^12    ^ 

E[XY]      =      -rx —. r4^7 r-^ r  + 


(X^-KX^2)(A2  +  X^2)  (A^+X^2^(X2  +  A^2^ 


and  the  correlation 


X,    -(X  +X   )6 

corr(X,Y)   =   p     =  -f^  e   -^   -^  (II-C-4a; 

x,y       X 


where 


X   =   X^  +  X2  +  X^2 
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If  6  =  0,  then  the  equation  (II-C-4a)  reduces  to 


p     =  -4^  (II-C-4b) 

xy       A 


This  characterization  can  be  represented  in  terms  of 
independent  random  variables.   For  6    >_  0 ,    if  there  exist 
independent   exponential  random  variables  U, ,  U^ ,  U_  and  U. 
with  respective  parameters  X,,  ^2'    ^lo'    '^12'  then 


X   =   min(U  ^U^) 


Y   = 


min(U  ,U  -6)     if    ^^  6 


L  min(U2,U^)       if   U^  <  6 


It  can  be  verified  formally  from  the  relation 

P[X>x,Y>y]       =      P  [U^  >  x,U2  >  y]  {P  [U^  >  x,U^  >  y+5  |U     >  5] 
X  P[U^  >  6]  }    +    P[U^  >  y]P[U     >  xlu     <  5]P[U2  >  6]  } 

In  the  5=0  case,  the  representation  yields 


X   =   min(U,,U  ) 


Y   =  min(U2,U2) 
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This  bivariate  distribution  has  a  line-discontinuity,  if 
E[X]  =  E[Y],  along  the  line  x  =  y,  i.e.,  in  the  case  where 
X  =  y  =  U-..   Thus  this  bivariate  exponential  is  not  as  smooth 
as  others.   In  particular  the  NEAR(l)  process  of  Lawrance 
and  Lewis  generates  bivariate  exponentials  with  a  continuous 
density  function.   But  if  we  consider  the  case  6  >_  0 ,  the 
line-discontinuity  will  be  removed.   In  either  case,  the 
resulting  (X,Y)  pairs  always  have  positive  correlation  as 
shown  in  equations  (II-C-4a)  and  (II-C-4b) .   If  we  use  these 
methods  to  generate  positively  correlated  random  vectors 
(X,Y)  for  exponential  marginal  distribution  with  unit  mean 
and  given  p,  we  have  to  compute  X,,  A^  and  X-«  for  given 
correlations  to  generate  independent  exponential  random 
variables.   For  simplicity  we  will  consider  the  method  with 
6  =  0.   From  the  relationship  '^  =  -^l  ■•"  ^^2  "*"  '^l  2' 


E[x]  =  Y-^j—    =  1 


E[Y]   =     /      =   1 
^2   ^12 


X 


P 


12 


X   ' 


we  can  compute  X,,  X^  and  X,^  as  functions  of  p 


^12   ^   2p/(l +  p) 


^1   ~   ^2   "   ■"■  "  ^12  • 
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Generating  Procedure 

1.  (Initialization) .   Compute  A  ,  A  and  X   for  given 
0  <_  p  <^  1. 

2.  Generate  three  independent  exponential  random  variables 
E,  ,  E^   and  E-.  with  unit  mean, 

3.  Rescale  the  independent  random  variables  and  set  as 


"i  "  ^i/*i 


^2   *   V^2 


U3   -  £37x^2 


4.   Define  X  and  Y  as 


X  =   min[Uj^,U  ] 


Y   =   min[U2,U2] 


5.   Deliver  (X,Y)  and  go  to  2  until  a  sufficient  number 
of  bivariate  pairs  have  been  generated. 
The  program  is  listed  in  the  appendix  and  the  scatter  plot 
and  bivariate  histogram  of  resulting  random  vectors  are 
shown  in  Section  V-D. 


38 


4 .   Gaver ' s  Negatively  Correlated  Bivariate  Exponential 
This  negatively  correlated  bivariate  exponential 
is  generated  from  the  following  considerations.   Suppose  a 
particular  system  has  N  defective  elements;  here  we  assume 
N  has  geometric  distribution  with  probability  density 
function 


f(x)   =   p(l-p)  ~  ,   X  =  1,  2,  ...  . 


The  time  to  failure  of  each  defective  element  is  T .  with 

1 

density  function  F(t).   Further  suppose  repair  is  carried 
out  after  all  defective  elements  are  discovered.   If  we 
define  R  as  repairing  time  which  is  distributed  as  G  (x) , 
provided  N  =  n,  and  define  T  as  minimum  of  T. /  then  the 
joint  distribution  function  of  T  and  R  is  obtainable  from 


P[T>t,R^x]       =         I    (1  -  F(t))"   P      G^(x) 

n=l  " 


where 


p^      =      (1-p)    p^"    ,      0    <    p    <    1,      n  =    1,    2, 


Then 


E[R|T>t]      =  ^^^^ 


1    -   P[l  -  F(t)  ] 
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where   E[R]    is   the    expected  value   of   an   element  repair   time 
Moreover, 


prRlT-M       -      FfRl    1   +   P[l-  F(t)  ] 
E[R|T-t]       -      E[R]     1    _    p[l-F(t)] 


=     E[R]    [1  +  i4rp(l-  <}>(t))] 


where 


<(>(t)      =      P[T<t]      =  ^^^^ 


1   -    P[l-  F(t)  ] 


From  this  regression,  we  know  that  if  T  is  long,  R  is  short 
and  vice  versa;  the  negative  relationship  is  clearly  present. 
Further  we  can  calculate  Cov[R,T)  as 


Cov[R,T)   =   P  E[R]  E[T](1  -  ^^^j^y^)   <   0, 


where  T(l)  has  the  distribution  of  the  smallest  of  a  sample 
of  two  from  cj)  (t)  . 

It  will  not  be  shown  that  one  can  select  F(t)  =  l-F(t) 
in  such  a  way  as  to  induce  exponentially  distributed  time 
to  system  failure,  T.   Since  we  know  that  exponential  G, 
geometrically  compounded,  yields  exponential  R,  the  outcome 
will  be  a  (T,R)  pair  with  exponential  marginals  and  negative 
correlation. 
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If 


F(t)   =   1  -  F(t) 


satisfies 


I       [F(t)]''  p""  ^    (1-p)   =   e 
n=l 


then  T  is  exponential  with  mean  1/A.   Solution  for  F 
yields. 


^<t'   =     ^     ,,     \        At       ^i° 

p  +  (1  -  p)  e 


which  is  a  logistic  distribution,  left  truncated  at  t  =  0 
If  X  =  1,  then  T  is  exponential  with  unit  mean,  while 
^^1{1)  ^    =  1/2,  so 


corr(R,T)   =   -  ^  E [R]  E [T]  . 


If  G  is  chosen  so  as  to  make  R  exponential  then  it  may  be 
shown  that 


corr(R,T)   =   -  |  =   -  i[l  -  ^1 


Consequently,  a  greatest  lower  bound  for  the  correlation  in 
this  model  is  -  y. 
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Generating  Procedure 

1.  If  0  <  p  <  -0.5,   Set  P  ^  -2p 

2.  Generate  a  geometric  random  variable  N  with  parameter  1-P 

3.  Generate  a  uniform  (0,1)  random  variable  U  and  define 
Y  as 


4.  Generate  a  gamma  random  variable,  Z,  with  shape 
parameter  n,  and  mean  is  n(l-p) 

5.  Deliver  (Y,Z)  and  go  to  2  until  a  sufficient  number 
of  pairs  are  obtained. 

For  obtaining  negatively  correlated  bivariate  exponential 
random  vectors,  this  algorithm  is  very  simple  and  one  of  the 
few  available.   Moreover  the  correlation  is  known  explicitly. 
The  program  is  listed  in  the  appendix  and  scatter  plot  and 
bivariate  histogram  of  resulting  random  vectors  are  shown 
in  Section  V-D. 

5 .   Arnold's  Vibariate  Gamma  Generator 

Arnold  (19  67)  developed  a  Trivariate  reduction  method 
to  generate  bivariate  gamma  random  vectors  having  positive 
correlation  with 


1/2 
corr  p  <  min (a, , a^) / (a,  a^) 


42 


where  a,  and  a-  are  the  given  shape  parameters  for  the  marginal 
distributions.   Letting  "gamma  (a, 3)"  denote  the  gamma  dis- 
tribution with  shape  parameter  a  and  scale  parameter  3/ 
so  that 


^<^'   =   BtT^  <f'""^  eKp(-x/B), 


where  r(a)  is  the  gamma  function  and 


2 

E[x]       =       a3         .  VAR[x]       =    a3     / 


the  trivariate  reduction  algorithm  proceeds  as  follows  for 
given  shape  parameter  values  a, ,  a^  and  p  to  generate  unit 
scale  gamma  variate  i.e.,  3  =  1. 

Generating  Procedure 

1/2 

1.  Generate  G,  ~  gamma  (a,  -p(a,a2)    /I) 

2.  Generate  G^  ~  gamma  (a2-p(a,a2)    ,1) 

1/2 

3.  Generate  G^  ~  gamma  (p(a,a2)    ,1) 

4.  Define  Y  and  Z 


y  <-     G^  +  G3 


z  -  G2  +  G3 
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These  Y  and  Z  are  unit  gamma  variates  with  6=1,  but  Y  and 
Z  can  be  multiplied  by  3-1  and  &^,    respectively,  to  obtain 
any  desired  scaling.   In  the  algorithm,  G-.  is  a  common  com- 
ponent of  both  Y  and  Z,  thus  inducing  ponitive  correlation. 

When  a  positive,  but  not  extreme,  correlation  is 
needed,  trivariate  reduction  is  an  excellent  algorithm, 
using  univariate  gamma  generators.   If  the  marginal  distri- 
butions have  the  same  shape  parameter  value  a,  then  the 
correlation  limitation  is  0  <_  p  <  1.   But  in  the  other  case 
the  upper  limit  is  smaller  than  1. 

Note  that  this  algorithm  exploits  the  infinite 
divisibility  of  the  marginal  gamma  variables  and  is  applicable 
to  any  infinitely  divisible  marginal.   It  is  commonly  used, 
for  instance,  to  get  bivariate  Poisson  random  variables. 
6 .   Schmeiser's  Bivariate  Gamma  Generator 

Schmeiser  (19  79)  developed  a  family  of  algorithms, 
any  member  of  which  can  generate  bivariate  gamma  random 
vectors  having  any  shape  parameters  a, ,  a_  and  allowable 
correlation  p.   In  this  section  we  will  discuss  his  general 
procedure. 

If  Z,  W,  and  VI      are  independent  gamma  random  varia- 
bles with  shape  parameters  r,  6,,  and  6^,  respectively,  and 
U  is  an  independent  uniform  [0,1]  random  variable,  and  either 
V  =  U  or  V  =  1-U;  then 


X,   =   F~-'-(U)  +  Z  +  W,  (II-C-5) 

J.      n  1 
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x^  =  f"-'-(v)  +  z  +  w^ 


where  F   and  F   are  the  inverse  distribution  functions 
n       m 

of  gamma  (n,l)  and  gamma  (m^l),  respectively,  distribution 
function.   The  resulting  X  ,  x^  are  gamma  random  variables, 
with  shape  parameters 


a.   =   n  +  r  +  6  ^  , 


a-   =  m  +  r  +  6  , 


(II-C-6) 


respectively.   This  result  follows  immediately  from  the 
reproducibility  property  of  the  gamma  distribution  and  from 
noting  that  F   (u)  and  F   (1-u)  are  each  random  variables 
having  CDF  F.   This  scheme  generalizes  the  trivariate 
scheme  but  brings  in  the  inverse  CDF. 

The  correlation  coefficient  is  defined  as 

E[f"-'-(u)  f'-'-Cv)  ]  -  nm  +  r 
_     n m /TT_p_7^ 

1   2  Ka^a^) 

The  remaining  problem  is  to  select  values  of  the  five 
parameters  n,  m,  r,  6,  and  ^^  to  obtain  the  desired  marginal 
distributions  and  correlation.   The  following  conditions 
must  be  satisfied. 


n  +  r  +  6^   =   a. 
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m  +  r  +  5    =   a 


E[F~^(u)  f"^(v)  ]  -  (nm)  +  r 

^      ^         =   p        (II-C-8: 


/     ^1/2 


H^,  H^/  r,     6^,     ^2      -      ^ 


Since  we  are  using  five  variables  to  satisfy  three  equality 
conditions,  finding  a  set  of  parameter  values  corresponds 
to  finding  a  feasible  solution,  rather  than  an  optimal  solu- 
tion, to  a  nonlinear  programming  problem. 

An  efficient  solution  procedure  for  determining 
parameter  values  is  important,  since  substantial  computation 
is  required  to  determine  whether  or  not  conditions  (II-C-8) 
are  satisfied  for  given  parameter  values.   Most  of  the 
computation  is  involved  in  calculating 

E[f"^(u)  f"-^(v)]  -  nm  -  r 

_     n m 

P   "         ,     .1/2 
(a^a2) 

since  the  expected  value  must  be  calculated  numerically 
using  any  one  of  the  following  three  integrals: 


^  -1      -1 
J  F^-^(u)  F^-^(u)  du  (II-C-9a) 


I  fJ^^(F^(x))  x^  exp(-x)  dx/r(n)         (II-C-9b) 
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/  F^-^(F^(-ln  y))(-ln  y)  ^  dy/r  (n)        (II-C-9c) 


if  p  >  0.   If  p  <  0,  then  replace  f"  (u)  with  F   (1-u)  in 

f      n  n 

equation  (II-C-9a) /  replace  F  (x)  with  1  -  F  (x)  in  equation 

(II-C-9b) ,  and  replace  F  (-In  y)  with  1  -  F  (-In  y)  in  equation 

n  n 

(II-C-9C) . 

Schmeiser  seemed  to  have  best  results  in  terms  of 
a  subjective  trade  off  between  speed  and  accuracy,  using  a 
24  point  Gaussian  method  for  the  integration.   He  selected 
the  parameter  values  from  the  feasible  values  satisfying 
conditions  (II-C-8)  by  making  the  curves  of  regression 
E[X, 1X2]  and  E[X2|X, ]  behave  as  desired. 

Generating  Procedure 

1.  Generate  X,  ~  gamma  (n,l) 

2.  U  ^  F  (X,  ) 

n   1 

Ifp  <  0,  U^l-u 

3.  Generate  Z  ~  gamma  (r,l) 


4.   Generate  W   -  gamma  (6,,1) 


5.   Generate  W^  ~  gamma  (6^,1) 


6.   X,  ^  X^  +  Z  +  W 


^2  ^  ^h"^^^^  "^  ^  ■*■  ^2 
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Based  on  these  procedures  he  developed  a  family  of  algorithms 
which  can  provide  variates  having  any  theoretically  possible 
correlation  p.   Anyway,  for  gamma  marginal  distributions, 
not  all  correlations  are  consistent  with  particular  shape 
parameter  values.   Schmeiser  shows  the  obtainable  correlations 
as  a  function  of  a^,    given  a,  =  1  and  5  as  shown  in  Section 
II-A.  Note  that  only  when  a,  =  a_  is  it  possible  to  obtain 
p  =  1-  Likewise  p  =  -1  is  not  possible  except  in  the  limit 
as  a,  and  a^  tend  to  infinity.  The  m.aximum.  and  minimum 
possible  correlations,  given  in  Moran  (1969),  occur  when 


^2  =  ^;J'\(^i" 


and 


X.   =   F  ^(1  -  F   (X, )) , 
2       a^  .      a^   1 


respectively,  where  F  (X)  and  F   (U)  are  the  cumulative 

a  a 

distribution  function  and  inverse  CDF,  respectively,  of  the 
gamma  distribution  with  shape  parameter  a. 
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III.   GENERAL  MIXTURE-TRUNCATION  METHOD 

Denote  by  (y,Z)  the  bivariate  random  pair,  where  each 
has  identical  marginal  continuous  distribution  F(x),  and 
denote  a  general  random  variable  from  this  distribution  by 
X.   The  argimient  is  not  specific  to  continuous  random  varia- 
bles; this  aspect  comes  in  only  in  the  computation  of  the 
correlations  and  can  be  developed  in  a  parallel  fashion  for 
discrete  marginal  distributions.   Let 


P   = 


1  -  a. 


l-a2    OL^ 


with  stationary  vector 


1- 


TT    =    TT  P    = 


1-  a. 


1-a,  +l-a^'  l-a,+l-a^ 


)   =   (tt-j^/  TT2)  , 


and  let  X-,  be  an  X  truncated  to  the  left  of  a  fixed  point  x  , 
X-  be  an  X  truncated  to  the  right  of  x  ,  so  that 


F^  (x)   =   P[X^£x]   = 


F(x) 
F(x^) 


if    X  <  x_ 
—  o 


if   X  >  X 
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F^  (x)   =   P[X  <x]   = 
^2 


r   0  if   X  <  X 

—   o 


F(x)  -  F(x^) 

O  ■  r- 

r =^7 r     if   X  >  X 

1  -  F(x^)  o 


In  addition  we  set  tt,  =  F(x  ),  tt-  =  1-  -it,  and  choose 

1      o     2       1 

Y  and  Z  as  follows. 

i)   Choose  Y  from  X-,  with  probability  tt,  and  then  choose 
Z  from  X,  with  probability  a,,  or  from  X^  with 
probability  1  -  a, . 
ii)   Choose  Y  from  X-  with  probability  tt^,  where  tt,  +  tt^  =  1, 
and  then  choose  Z  from  X,  with  probability  1-  a„,  or 
from  Xp  with  probability  a^. 

If  we  choose  (Y,Z)  as  in  the  above  procedure,  then  we  can 
make  the  following  two  theorems. 

A.   MARGINAL  DISTRIBUTIONS 
Theorem  1. 
The  marginal  distribution  of  (Y,Z)  becomes  F(x)  for  both 

Y  and  Z. 

Proof 

1.   Marginal  distribution  of  Y 

By  definition  Y  is  the  mixture  of  X,  and  X^  with  proba- 
bility 7T^,  TT^  respectively.   That  is 


Fy(x)   =   ^^  F^^  (x)  +  Tr2  F^  (x) 
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F  (x)   = 

y 


^1  ¥WJ  '  '2   '   ' 


F(x)  -  F(x  ) 
^1  •  ^  ■"  ^2    1-F(x J 


if   X  <  X 


if   X  >  X 


But  since  we  define  tt,  =  F(x  )  ;  tt_  =  1  -  -rr,  =  1-F(x  ),  we 

1      o    2       1         o 

have 


Fy(x) 


F(x  ) 

7T,   — -^    =    F(X) 

1    TT  - 


F(X)  -  TT. 


^1  ■"  ^2 


if   X  <  X 

—  o 


=   F(x)    if   X  >  X 


=   F (x)     in  all  cases. 


So  Y  has  the  marginal  distribution  F(x) 


2.   Marginal  distribution  of  Z 
If  Y  is  from  X,,  then 


F„  (x) 
^1 


a^  F^  (x)  +  (1  -  ai)F^  (x) 


^IHCT^  ^^-^1^  •  ' 


if  X  <  X 
—  o 


F(x)  -  F(x^) 
a^  .  1  +  (1-a^)   1  ,  p(^  )    if  ^  >  ^o 
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If  Y  is  from  X  ,  then 


F„  (X) 
^2 


(1  -  a^)  F^  (x)  +  a^  F^  (x) 


^^-^2^  ¥WT^    ^2    '    0 


F(x)  -  F(x^) 
^^-^2^  •  1  ■"  ^2    1-F(x  ) 


if   X  <  X 
—  o 


if   X  >  X 


So, 


F^CX)    =    TTj^  F^   (X)   +  TT^  F^   (x) 


.  .   F(x)   ,     ^     (.  .    F(x) 


F(x)-  FCXq) 


if  X  <  X 


L  TT^(a^+  (1-a^)     1  _  F(x  ) 


if  X  >  X 


F(x)  -  F(x^) 
+  Tr2((l-a2)  +  a^    i-f(x  ) 


and  we  defined  tt,  and  tt^  as  follows. 


1  - 


1  -  a,  +  1  -  a^ 


1  -  a. 


l-a^+l-a^ 
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From  this,  we  know 


1  -  a    =   -i  (1  -  a,) 

2        TT  1 


If  we  use  this  relationship,  then  F_(x)  =  F(x)  in  both 
cases.   So  Z  also  has  the  marginal  distribution  F(x).   The 
result  is  a  consequence  of  the  fact  that  _n_  is  defined  to  be 
the  stationary  vector  associated  with  P. 

B.   THE  PRO DUCT -MOMENT  CORRELATION 
Theorem  2. 
The  correlation  coefficient  between  Y  and  Z  becomes 

P   =   3  M  , 


where 


and 


where 


-1  ^  6  =  a^^  -  (1  -  a2)  1  1   , 


M   =  i^Sj^-M^)       Tr^TT2/a^ 


r  °  V  H  g'(x) 
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F(x)  -    F(x    ) 

^2      =  ^    ^   *^    1   -  F(x    )° 

X  o 

o 


X 

2  r    °      2     .    F(x)  2 


^2        =         J         ^      ^      1    -    F(x    ) ^2 


a^^      =         /    x^    dF(x)    -    E[X] 
^  0 


=      0^2    TT^   +    02^    TT^    +    (y^-  ^2)^    ^1    ^2 


Proof 


cov[Y,Z] 


E  [YZ]    -    E[Y]    E[Z] 
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E[YZ]       =      E     [E[YZ  |Y,Z    £    S]  ] 


where 


Further, 


=      E[YZ|Y    e    X^,Z    e    X^]P[Y    z    X^,Z    e    X^] 
+    E[YZ|Y    e    X^,Z    e    X^]F[Y  e  X^,Z    e    X^] 


+   ELYZlYeX^/Z    e    Xj^]P[Y    e    X^,Z    e    X^^] 


+    E[YZ|Y    e    X^/Z    e    X^]P[Y    e    X^ / Z    e    X^] 


=       E[X^]     E[X^]     TTj_aj_    +    E[X^]     ElX^]     TT^(l-a^) 


+    E[X2]    E[X^]     7T2(l-a2)     +    E[X2]    ELX^]     i^^o.^ 

2 

:       y,        a^    TT^    +    y^    y2(l-a^)     tt^^ 

+    y^    y2(l-a2)Tr2    +    y2       a2Tr2 


yj_      =       E[X^],  y2       =       E[X2] 


E[Y]       =      Eg[E[Y|Y    £    S]] 


E[Y|Y    £    Xj^]P[Y    £    X^]    +    E[Y|Y    £    X2]P[Y    z    X^\ 


E[X^]     TT^    +    E[X2]     TT2 


^1    ^1    "*"    ^2    ^2      ~      E[Z] 


by  Theorem  1. 
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Also  by  Theorem  1, 


o^      =   E[Y^]  -  (E[Y])^ 


2        2 


If  we  put  together  these  formulae  into 


E[YZ]  -  E[Y]E[Z] 
^Y,Z  a^  o^ 


we  get 


2 

(y^  -  y^)  TTj^TT  (a^  -  (l-a2)) 

Py,Z   =   2 

^X 


and  let 


a,  -  (1  -  0.2) 


M  =   (y^  -^2^  ^1^2/^X^ 


then 


p   =   6  M  . 

C.   GENERAL  ALGORITHM 

We  give  here  three  algorithms  for  implementing  the 
bivariate  mixture-truncation  method,  which  we  call  the  FXO 
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method,  the  UXO  method,  and  the  TXO  method.   All  of  these 
methods  are  exactly  the  same  except  in  how  the  algorithm 

chooses  X  ,  the  truncation  point,  from  the  x   range  [x  ,x  ] . 

o  c  '  o  2,   u 

The  first  procedure,  called  the  FXO,  choose  x  as  a  fixed 
point  of  X  and  x  and  uses  the  same  x  during  the  entire 
routine.   The  second  procedure,  called  the  UXO,  chooses  x 
uniformly  from  [x  ,x  ]  and  repeats  this  step  in  every  call 
to  the  algorithm.   The  third  procedure,  called  the  TXO,  is 
the  same  as  the  UXO  procedure  except  in  that  it  uses  a  tri- 
angular distribution  instead  of  uniform.   It  is  necessary  to 

fix  these  choices  of  x_  becuase  in  general  there  is  more  than 

o  ^ 

one  X  which  will  give  a  bivariate  pair  (Y,Z)  with  the  given 

marginal  distribution  and  given  correlation.   The  first 

procedure,  FXO,  is  defective  in  terms  of  their  discontinuity 

of  distribution  while  the  second  and  the  third,  UXO  and  TXO, 

are  satisfactory  in  this  respect.   The  choice  of  the  midpoint 

of  the  interval  [x^,x  ]  for  x   in  FXO  is  based  on  experience 

£.   u       o  '^ 

presented  later  for  various  marginal  distributions.   Note 
that  the  algorithm  described  here  is  inefficient  in  that  it 
generates  the  truncated  variables  X,  and  X^  by  comparing  ran- 
dom variables  X  to  x  until  one  which  is  respectively  greater 

than  or  less  than  x   is  found.   More  efficient  methods  can 

o 

be  found  in  special  cases  such  as  the  exponential,  but  the 
present  algorithm  requires  only  a  generation  of  univariate 
random  variables  X  without  regard  to  the  method  used  to  do 
this.   Of  course  initialization  is  required  and  this  is 
specific  to  each  marginal  distribution. 
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General  Mixture-Truncation  Method 

1.  (Initialization) 

i)   For  given  marginal  distribution  F(x)  and  correla- 
tion coefficient  p  find  x   ranges  [x,,x  ] 

O      ^       £   U 

2.  Define  truncation  point  x 

^      o 

*  FXO  method 

i)   X    =   i(x,  +  X  ) 
o      2   JJ,    u 

*  UXO  method 

i)   Generate  a  uniform  [0,1]  random  variable  V^ 

ii)   X   =   x„  +  (x  -  X,  )  *  V, 
o       £.      u    Jl.      1 

*  TXO  method 

i)   Generate  two  uniform  [0,1]  random  variables 

ii)    X^   =   x^  +  Xj_  +  Xj 

where 

X    =  -^{x      +    X  ) 
m     2   £    u 


^1  =   (^m-^>  *  ^1 


x^   =   (x  -  x^)  *  V^ 
2       u   m     2 


3.   Compute  parameters  value,  tt  ,  t\    ,    a,,  a^ 


4 .   Choose  type  for  Y 

i)   Generate  a  uniform  [0,1]  random  variable  U 
ii)   If  U  <_  7T,  ,  go  to  9 
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5.   Y  is  an  X 

i)   Generate  a  random  variable  X  from  F(x) 


ii)   If  X  >  X  ,  set  Y  -«-  X  and  go  to  6 


iii)   Otherwise  return  to  5.i) 


6.   Choose  type  for  Z 


i)   Set  U  ^  (  (U-  TT^)/(1  -  7Tj_)  ) 
ii)   If  U  £  1  -  a2,  go  to  8 


7.  Z  is  an  X^ 

i)   Generate  a  random  variable  X  from  F(x) 
ii)   If  X  >  X  ,  set  Z  ^  X  and  go  to  11 
iii)   Otherwise  return  to  7.i) 

8.  Z  is  on  X, 

i)   Generate  a  random  variable  X  from  F(x) 
ii)   If  X  <_  X  ,    set  Z  ^  X  and  go  to  11 
iii)   Otherwise  return  to  8.i) 

9.  Y  is  an  X, 

i)   Generate  a  random  variable  X  from  F(x) 


ii)   If  X  <^  X  ,  set  Y  ^  X  and  go  to  10 
iii)   Otherwise  return  to  9.i) 


10.   Choose  type  for  Z 
i)   Set  U  -r-   u/tt 
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ii)       If    U    <_   a-  ,    go   to    8 
iii)      Otheinvise   go   to    7 


11.   Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  method  until  a  sufficient 
number  of  random  vectors  are  obtained. 

D.   BIVARIATE  DISTRIBUTION  FUNCTIONS 

From  Theorem  1,  in  Section  III-A,  we  know  that  if  Y 
is  from  X,,  then 


F2(z|y) 


a^  F^  (z)  +  (1-  a^^)  F^  (z) 


F(z) 
■1  F(x^) 


F(z)  -  F(x^) 
^'l'-  (l-«l)    1  -  F(x  ) 


if   z  <  X 


if   z  >  X 


and  if  Y  is  from  X   then 


F^iz\Y) 


(1  -  a^)     F^    (Z)  +  a^    Fx  (z) 


(1-  aj 


F(2) 
2'    F(x^) 


F(z)  -  F(x  ) 
vx  ^^j  "^2   1  -  F(x  ) 


if   z  <  X 
—  o 


if   z  >  X 
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By  using  these  we  can  define  the  bivariate  distribution 
function  as  follows. 

F(y,2)   =  V[Y    <_Y,Z    <_  z] 

j  P[Z  <_  z  |y  =  u]  dP[Y  j^  u] 


where 


P[y  <  u]   =   F(u) 


P[Z  <  z  |y  =  u] 


F(z)  -  F(x^) 
^1^   (l-«l)    1-F(x  ) 


(1-  aj 


F(z) 
2'    F(x^) 


F(z)  -  F(x^) 
^^-'^2^    -^  ^^2   1-F(x  ) 


if  u  <  X  ,  z  <  X 
—  o'   —  o 


if  u<x  /  z>x_ 
—  o      o 


if  u  >  X  ,  z  <  X 
o    —  o 


if  u  >  X  ,  z  >  X 
o      o 


So,  if  we  put  these  together,  integrating  with  respect  to 
dP  [Y  <_  u]  ,  we  get  the  final  result: 
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P[Y<y,Z<z]     = 


F(z)F(y) 


if  y  <  X  ,   z  <  X       (III-D-1) 
^  —  o       —  o 


(1-a   ) 

{a,  + ^[F(z)-F(x  )]}F(y)    if  y  <  x  ,   z>x        (III-D-2) 

Itt-  O  —    O  O 


(l-aj 
{a,  + ^[F(y)-F(x J]}F(z)    if  y  >  x  ,   z  <  x^      (III-D-3) 

J.  7T~  O  O  —     O 


Vl  ^   ^1-^1^  [F(z)-F(x^)]   ^ 


(III-I>4) 


+{(l-a2)  +^[F(z)-F{xQ)]}[F(y)-F(x^)]    if  y  >  x^,   z  >  x^ 


For  example,  the  expression  (III-D-4)  is  obtained  as 


F(y,z)       =  /    P[Z    <_  z|y   =   u]    dF(u)         y    >    x    ,       z    >    x 

—  CO 


/      P[Z^z|y   =    u^x]    dF(u! 


+         /P[Z<_z|Y   =    u>x]dF(u) 

X 


(1-  a    ) 
{a^    +   — ^[F(z)    -    F(x^)]}    F(x^) 

±  7T  ,,  O  O 


+       (1-a^)    +   — [F(z)    -    F(x^)]}[F(y)    -    F(x^)] 

2  TTo  O  O 


It   is   easily   seen   from    (III-D-2)    that  when   z   ^   «>, 
P[Y    <_  y,Z    <_  oo]    =   F(y);    from    (III-D-3)    that   as   y   ->    «., 
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P  [Y  <_  °°,  Z  <_  z]  =  F(z)  and  from  (III-D-4)  that  as  y  ^  <=° , 

P[Y  <_  ~,Z  <_  z]    =  F(z)  and  that  z  ^  «>,  P[Y  <_  y,Z  <_  oo]  =  F(y)  . 

In  particular  from  (III-D-4)  we  have  that,  as  y  ^  °°, 


F(y,z)   ->   ct^F(x^)  +  (1  -  a2)  [F(z)  -  F(x^)  ]  +  (l-a2)Tr2 


+  a^lFCz)  -  F{x^)] 


=   ct  F(x  )  +  (1-a  )tt   +  [F(z)-F(Xo)] 


•1"  '^O 


2'  2 


F(z)  +  (l-a2)TT2  -  (l-a^)TT^   =   F(z) 


where  at  the  last  step  we  used  the  facts  that  F(x  )  =  tt, 
and  Tr,(l-a,)  =  TT^d-ct^).   If  F(x)  is  absolutely  continuous 
with  probability  density  function  f(x)  then  the  joint  p.d.f. 
for  the  bivariate  pair  (Y,Z)  is 


f(y,2)   = 


a,  (1  -  a, ) 
TT2(l-a2) 


1-a, 
1-a. 


f(y)  f(z) 


f(y)  f(z) 


^2 


if  y  <  X  ,  z  <  X 
■^  —  o    —  o 


if  y  <  X  ,  z  >  X 
•^  —  o       o 


if  y  >  x^,  z  £  Xq 


if  y  >  X  ,  z  >  X 
^  o       o 


(III-D-5) 


(III-D-6) 


(III-D-7) 


(III-D-8) 
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Note  that  there  is  a  discontinuity  in  the  density  function 
as  one  crosses  the  boundaries  of  the  four  quadrants  defined 
by  the  lines  y  =  x  ;  z  =  x  .   The  density  is  the  same  in  the 
first  and  third  quadrants.   The  multipliers  of  f(y)f(z)/7T- 
are  the  same  in  all  four  quadrants  iff  there  is  independence 
This  occurs  when  (1-a,)  =  a^  so  that  it,  =  cto  ^^*^ 
f(y,z)  =  f(y)f(z)  for  the  whole  range  of  y  and  z. 
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IV.   BIVARIATE  UNIFORM  GENERATOR 

The  uniform  random  variable  is  a  continuous  random 
variable  with  probability  density  function  which  is  constant 
over  the  interval  (a,b)  and  zero  otherwise;  the  density 
function 


f  (X) 


b-a' 


if   a  <  X  <  b  , 


0;       otherwise  ; 


mean  =   E[X]   = 


b+a 


2  ' 

2 


variance   =   VAR[X]   =  S^  ~  ^^ 


12 


We  note  that  if  U  has  a  uniform  distribution  over  [0,1] 
then  X  =  a  +  (b-a)U  has  a  uniform  distribution  over  [a,b] . 
So  we  will  only  be  concerned  with  algorithms  for  generating 
uniform  [0,1]  distributions. 

Two  algorithms  for  generating  uniform  bivariate  pairs 
are  given  in  Section  II,  Gaver ' s  and  Lawrance  and  Lewis's, 
and  since  they  have  relatively  smooth  distributed  functions  and 
are  simple  to  generate,  they  are  probably  preferable  to  the 
uniform  bivariate  random  variable  generated  by  the  mixture- 
truncation  method  unless  x   is  taken  to  be  random  (of  course 
the  correlation  for  the  Gaver  bivariate  uniform  is  difficult 
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to  compute) .   However  the  development  in  this  section  of 
the  uniform  mixture-truncation  method  does  help  to  fix  ideas 
on  the  problem  of  initialization  and  determination  of  the 
possible  range  of  correlations  attainable  in  other  cases. 
To  use  the  mixture- truncation  method  for  generating  bivariate 
random  vector  (Y,Z)  whose  marginal  distribution  is  uniform 
over  [0,1]  identically  and  has  given  correlation  q,    we  should 
use  uniform  [0,1]  distribution  functions  as  F(x).   Then  X,, 
which  is  X  constrained  to  be  on  the  left  side  of  the  trunca- 
tion point  X  ,  becomes  uniform  over  [0,x  ]  and  X^ ,  which  is 
constrained  to  be  on  the  right  side  of  x  becomes  uniform 
over  [x  , 1] . 

A.   DETERMINATION  OF  PARAMETERS  IN  THE  UNIFORM  MIXTURE- 
TRUNCATION  METHOD 

Because  X^  is  uniform  [0,x  ]  and  x_  is  uniform  [x  ,1] / 


E[X^]   =   p^   =   ^ 


2 
VAR[X^]   =   a^^   =   -^ 

1  +  X 
E[X2]   =   ^2   =  -^ 


2      ^l-^o^^ 
VAR[X2]   =   a^"   =   j2^- 


and 


1  -  a^ 

TT^   =   F(x^)   =  -^ -^ =   x^  ,  (IV-A-la) 

1         O       l-a,+l-a,,       O 


66 


1  -  a 
TT,   =   1  -  F(x^)   =   .j -4 =   1  -  x^   (IV-A-lb) 


If  we  use  these  formulas  in  Theorem  2  of  Section  III,  we  get 


M      =      3   X    (1  -  X   ) 
o  o 


°1   -   ^o 
1    -    =^0 


p       =       3    X    (a,     -    X    )  (IV-A-2) 

o      1  o 


Then   from    (IV-A-2) ,     (IV- A- la)    and    (IV-A-lb) 


a,       =      .5^-  +    X  (IV-A-3a) 

1  3x  o 


""l 
a^      =       1 ^(1  -  a.) 

2  Tr2  1 

2  ^ 
X         +  X 

=      1 +   ^ (IV- A- 3b; 


1-x^  3(1  -x^) 

o  o 


Furthermore,  a,  and  a^  are  probabilities  so  they  have  to 
satisfy  the  following  conditions: 


0  <_  a^  <_  1, 


0  1  ^2  -  ■'■  * 
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By  solving  these  two  inequality  equations,  we  can  get  the 

feasible  range  of  x   for  given  p  as  follows. 

o 


x„      =      Lower  bound  of   x 
Jii  o 


1/2    -    /I/ 4  -(1/3)  p  if      p    >    0 


/-p73 


(IV-A-4) 


if      p    <    0 


X        =      Upper  bound  of   x 
u  ^^  o 


1/2    +    va/4  -a/3)  p  if      p    >    0 


1   -    /-p73 


(IV-A-5) 


if      p    <    0 


After  finding  x„  and  x  ,  choose  x  within  the  range  [x, ,x   ] 

I  u  o  ^     £   u 

by  way  of  either  the  FXO,  UXO,  or  TXO  method.   And  by  formulas 

,(IV-A-la),  (IV-A-lb) ,  (IV-A-3a)  and  (IV-A-3b) ,  we  can  get 

TT,,  TT^f  a,  and  a^.   Also  we  can  expect  some  limitations  in 

p  because  of  x  .   From 
o 


M  =   3  x^(l-  x^) 
o     o 


-1  <_  3  =  a,  -  (1  -  o.^)    1  1 


P 


M 
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We  can  compute  that 


p      =   3/4   if  Q    =    1,      X  <  X   =  1/2  <  X 
max      '  '   £  —  o    '   —  u 


p  .    =   -3/4  if   3  =  -1,  x  <  X  =  1/2  <  X 
mm       ^  i  —     o  '   -  u 


In  another  way,  note  that  if  one  chooses  x  =  1/2  and 

a^  =  ot-  =  1,  from  (IV-A-2)  we  attain  the  highest  allowable 

correlation  p  =  3/4,  and  when  a^  =  0,  p  =  -3/4,  the  lowest 

attainable  correlation. 

It  turns  out  that  for  this  choice  of  x  =  1/2  (in  the  FXO 

o    ^ 

method) ,  the  bivariate  uniform  distribution  is  also  about 
the  smoothest,  as  can  be  seen  in  a  later  section  (IV-D), 
and  from  the  bivariate  distribution. 


B.   GENERATING  PROCEDURE 

j    We  developed  here  all  of  the  three  procedures,  the  FXO 
method,  the  UXO  method,  the  TXO  method  for  generating 
bivariate  random  vectors  whose  marginal  distributions  are 
uniform  over  [0,1]  and  correlation  coefficient  is  p .   In 
this  algorithm  we  used  direct  generating  procedures  for  X, 
and  X^  instead  of  comparing  random  variables  X  to  x  until 
one  which  is  respectively  greater  than  or  less  than  x   is 
found. 
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Uniform  Mixture-Truncation  Method 
1.   (Initialization) 

i)   For  given  -3/4  <  p  <  3/4,  find  x^  and  x 


1/2  -  /1/4  -  1/3  p   if  p  >  0 


/^Jn 


if  p  <  0 


X 


u 


1/2  +  /1/4  -  1/3  p   if  p  >  0 


1  -  /-p/T 


if  p  <  0 


2.   Define  truncation  point  x  . 

^      o 


*   FXO  method 


1     X, 


(x^  +  x^)/2 


*  UXO  method 

i)   Generate  a  uniform  [0,1]  random  variable  U, 

ii)   X   =x+(x-x)*U, 
o       t  M.  I  1 

*  TXO  method 

i)      Generate   two   uniform    [0,1]    random  variables 

ii)       X        =      x„    +   X,    +   x^ 

o  i,  1  2 

where 


X,      =      (x^   -   X   )     *   V 

1  mil 
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x„      =       (x      -   x„)     *   V^ 
2  urn  2 


X        =       (x      +   xJ/2 
m  u  g^ 


3.      Compute    tt,,    tt^z    a,    and   a«    as 


^1      =      F(^o^      =      ^o 


^2      =      ^    "    ^1 


a,       =  -^  X      +    X 

1  3      o  o 


a,      =       1    -    -i(l    -    a.) 
2  7T2  1 


4.   Choose  type  for  Y 

i)   Generate  a  uniform  [0,1]  random  variable  U 


ii)   If  U  <_  TT,  go  to  9 


5.   Y  is  on  X2 

I      i)   Generate  a  uniform  [0,1]  random  variable  W, 

I 


ii)   Y  ^  X   +  (1.0  -  X  )  *  W, 
o  o     1 


6.   Choose  type  for  Z 

i)   Set  U  ^  ((U  -  7:^)/(l  -  TT^)) 

ii)   If  U  <  1  -  a^,  go  to  8 
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7.  Z  is  an  X^ 

i)   Generate  a  uniform  [0/1]  random  variable  W^ 

ii)   Z  ^  X  +  (1.0  -  X  )  *  W^ 
o  o     2 

iii)   Go  to  11 

8.  Z  is  an  X 


1 
i)   Generate  a  uniform  [0,1]  random  variable  W 


2 

ii)   Z  ^  x^  *  W^ 

iii)   Go  to  11 

9.   Y  is  an  Xj^ 

i)   Generate  a  uniform  [0,1]  random  variable  W, 

ii)   Y  ^  X   *  W, 

10.  Clioose   type    for   Z 
I  i)      Set   U  ^   U/tt 

ii)       If   U    <_  a,    go   to    8 

iii)   Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4,  for  the  FXO  method,  or  go 
to  2  for  the  UXO  method  and  the  TXO  method  until  a 
sufficient  number  of  random  vectors  are  obtained. 

The  programs  are  listed  in  the  appendix  and  the  scatter 
)lot  and  bivariate  histogram  of  resulting  random  vectors  are 
shown  in  Section  IV-D. 
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C.   REGRESSION  OF  Z  ON  Y  FOR  GIVEN  p 

Schmeiser  (19  79)  in  developing  a  bivariate  gainina  method 

fixed  the  free  parameter  by  specifying  the  regression  of  Z 

and  Y  and  vice- versa.   We  do  this  now  for  the  uniform  case 

with  fixed  x   and  uniformly  distributed  x  .   First  in  fixed 
o  -^  o 

x  case  we  have,  for  y  <  x 
o  "^  —  o 

E[Z|Y  =  y,y^x^]   =   a^E[X-^]  +  (1  -  a^)  E[X2] 


=   1/2  -  -^   p 
6x 
o 


and  for  y  >  x 
^  o 


E[Z|Y  =  y,y  >  x^]   =   (1-^2)  E[X^]  +  a^    ^^^2^ 

1 


=   1/2  +  p 


6(l-x^) 


These  are  the  step  function  and  clearly  not  linear  in  y. 

If  X  has  uniform  distribution  over  [x„,x  ],  we  have  the 
o  I'    u 

following  results.   These  results  are  dependent  on  Y  values 
If  y  ±  ^nt    then  we  have 


X 

u 
E[Z|Y  =  y]   =     /   E[Z|Y  =  y,X=x^,y£x^]  f(x^)  dx^ 


X 

r  ^  /I  1    ^    1     A 

1    ( o- -  7 —  p)  dx 

^            2  6x  X  -  x„    o 

X  o     u   e 
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E[Z|Y  =  y]       = (^  X     -^x.-|ln-^) 

'         ^  x-x„2u2      5i6  X 


If  Xrt    <   y    <    X    ,    then  we   have 


^£   :i  ^   _  ^u 


Y 
E[Z|Y  =  y]       =  J      E[z|Y  =  y,    X=x    ,    y>x^]f(x^)    dx^ 

x„ 


X 

u 

+       {        E[z|Y  =  y,    X=x    ,    y<x    ]    f(x    )    dx 
-'  *•'         -^  '  o—    o  o  o 


/"    'I  ^    6(l-x    )'    f '-o'    d^o 


X 

o      1 

+      f         (^  -   ^)    f  (x    )    dx 
-'  2         6x0  o 

y  o 


1/1  1  Pi     i~y       u    K 

(^  X    -  TT  x^   -  -7-  In  — ^  Tj ) 

x-x2u2£       6  yl-x 


If  y   >    X    ,    then  we 


X 

u 


E[ZlY  =  y]       =  I         E[zlY  =  y,    X=x^,    Y  >  ^o^    ^^^o^    ^^o 


^£ 


X 

u 

dx 
.-  .  -  o 

X 

I 


/   "    '?^    6(1 -X   )'    ^<-o' 


T  T  1  1    -    X 

1/1  1  Pi  Uv 


x-x2u2!J.       6  1-x 
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From  the  results  we  can  see  that  if  we  use  uniform  distri- 
bution for  X  we  can  get  smoother  regression  functions  if 
X.J  _<  y  _f_  X   although  we  still  have  step  functions  for  the 

y  <  x„  and  y  >  X  cases. 

■^  —  £       —  u 

D.   SIMULATION  RESULTS 

We  will  show  here  the  scatter  plots  and  the  bivariate 
histograms  of  the  mixture-truncation  bivariate  uniform  (0,1) 
random  vectors  with  correlation  p  =  0.3  and  p  =  -0.3  by  the 
FXO,  UXO  and  TXO  methods  in  Figures  (IV-a) ,  (IV-b)  and 
(IV-c) ,  respectively. 

In  addition,  we  will  show  the  scatter  plot  and  bivariate 
histogram  of  sample  size  2000  from  Gaver ' s  transformation  and 
Lawrance-Lewis  method  in  Figures  (IV-e)   and   (IV-d), 
respectively.   The  Gaver  method  does  not  have  known  correla- 
tion but  the  value  of  p  to  give  p  in  the  exponential  case 
is  used. 

As  we  mentioned  earlier,  the  FXO  method  has  some  discon- 
tinuity at  the  truncation  point.   But  the  UXO  and  TXO  methods 
generate  relatively  smooth  continuous  distributions.   In  this 
respect,  we  say  that  the  FXO  is  defective  and  the  UXO  and 
the  TXO  are  relatively  satisfactory. 

The  computed  correlation  from  subroutine  BIVHST  (Bivariate 
histogram)  is  a  little  different  from  given  correlation.   But 
this  can  be  assumed  as  a  sampling  error.   To  check  this  error, 
we  simulated  10  times  with  given  correlation  p  =  0.1.   From 
this  simulation  we  get  the  following  results: 
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p   =   mean  of  computed  correlation   =   0.105 


VAR[p]   =   variance  of  computed  correlation  =   0.0004 


a[p]   =   standard  deviation  of  mean   =   0.0064 
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Figure    IV-al . 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  uniform  with 
p  =  0.3.   Here  x   is  fixed  at  the  midpoint 
between  the  lower  and  upper  bounds. 
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Figure  IV-a2.  Bivariate  historgram  for  sample  of  size  2000 
from  mixture- truncation  bivariate  uniform 
with  p  =  0.3.   Here  x   is  fixed  at  the 
midpoint  between  the  ?ower  and  the  upper 
bounds . 
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Figure  IV-a3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  uniform  with 
p  =  0.3.   Here  x   is  fixed  at  the  midpoint 
between  the  lower  and  the  upper  bounds. 
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Figure   IV-a4 . 


Bivariate   histogram   for   sample  of   size    2000 
from  mixture-truncation  bivariate   uniform 
with   p    =  -0.3.      Here   x      is    fixed  at   the 
midpoint   between   the    lower   and   the   upper 
bounds . 
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Figure  IV-bl 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  uniform  with 
p  =  0.3.   Here  x   is  taken  to  be  uniformly 
distributed  between  the  lower  and  the 
upper  bounds . 
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Figure    IV-b2 . 


Bivariate   histogram   for   sample  of   size    2000 
from   mixture-truncation  bivariate   uniform 
with  p  =    0.3.      Here   x      is    taken   to  be 
uniformly   distributea  between  the   lower   and 
the   upper   bounds . 
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Figure  IV-b3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  uniform  with 
p  =  -0.3.   Here  x   is  taken  to  be  uniformly 
distributed  between  the  lower  and  the  upper 
bounds 
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Figure   IV-b4. 


Bivariate  histogram   for    sample   of   size    2000 
from  mixture-truncation    bivariate   uniform 
with   p    =   -0.3.      Here   x      is    taken   to  be 
uniformly  distributed  Between   the   lower   and 
the   upper  bounds . 
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Figure   IV-cl 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  uniform  with 
p  =  0.3.   Here  x  has  triangular  distribution 
between  the  lower  and  the  upper  bounds . 
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Figure   IV-c2 


Bivariate   histogram   for   sample   of   size    2000 
from  mixture-truncation  bivariate   uniform 
with   p    =   0.3.      Here   x      has   triangular 
distribution  between   the   lower   and  the 
upper   permissible   bounds. 
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Figure   IV-c3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture- truncation  bivariate  uniform  with 
p  =  -0.3.   Here  x  has  triangular  distribu- 
tion between  the  lower  and  the  upper  bounds 
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Figure   IV-c4. 


Bivariate   histogram  for   sample   of   size 
2000    from  mixture-truncation  bivariate 
uniform  with   p    =   -0.3.      Here   x      has 
triangular   distribution  between   the   lower 
and  the   upper   bounds . 
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Scatter  plot  for  sample  of  size  2000  from 
Lawrance-Lewis ' s  method  with  p  =  0.3  and 


a    = 


(2+e) 


(4  + 


1) 


3+(2+e)p 

and   S    is    uniform  over    [ 


2p 

3-p' 


1]. 
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Bivariate   histogram   for   sample   of   size    2000 
from  Lawrance-Lewis' s  method  with   p    =    0.3 
and 
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(2  +  6) 


3+(2+3)p     '6 


(f+    1) 


and   6    is    uniform  over       [- 


2p 
3-p  ' 
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Figure  IV-el.   Scatter  plot  for  sample  of  size  2000 

from  Gaver ' s  transformation  with  p  =  0.6 
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Figure    IV-e2.      Bivariate  histogram  for   sample  of   size 

2000    from  Gaver ' s   transformation  with 
p  =    0.6. 
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V.   BIVARIATE  EXPONENTIAL  GENERATOR 

Another  probability  distribution  of  major  interest  in 
simulations  is  the  exponential.   The  cumulative  distribution 
function  and  probability  density  function  for  the  exponential 
are,  respectively. 


F(x)   =   1  -  e"^^     X  >  0 


and 


f(x)   =  X    e~^^       X  >  0 


The  expected  value  of  the  exponential  distribution  is 


E[X]   =   1/X 


and  the  variance  is 


VAR[X]   =   1/A^ 


The  problem  of  generating  exponential  deviates  reduces  to 
one  of  generating  "unit"  exponentials,  i.e.,  those  with 
A  =  1.0,  and  then  multiplying  the  result  by  whichever  1/X 
is  necessary  to  give  the  desired  distribution.   That  is,  if 
the  random  variable  E  has  the  exponential  (X  =  1)  distribution 
then. 
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X   =   E  *  i 


nas  the  exponential  (A  =  n)  distribution.   In  this  section, 

lie   considered  only  unit  exponentials  as  marginal  distributions 

for  bivariate  pairs. 

k.       DETERMINATION  OF  PARAMETERS  IN  THE  EXPONENTIAL  MIXTURE- 
TRUNCATION  METHOD 

Because  X,  is  an  exponential  (A  =  1)  truncated  to  the 

left  of  X   and  X2  is  also  exponential  (X  =  1)  truncated  to 

the  right  of  x  ,  we  have 


X 


1  -  —  x^; 

TT,    O 


x 

VAR[X,]   =   0^2   ^   ^/  °  x'  d  |{^  _  .^2 


^2     2 
=   1 2  ^o  ' 

^1 


F(x)  -  F(x  ) 

ElX^l   =   P2   =     /   '"^   l-F(K) 

X  o 

o 


=   1  +  x^; 
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"   2    F(x)  -  F(x  )      2 
VAR[X  1   =    /   X  d    X  .  F(^  )  -  ^'2 
X  o 


o 
=  1; 

md  from  definition, 

1  -  a2 

7T. 


1      1  -  a-,  +  1  -  a. 


=  F(x^) 


-X 

=   1  -  e   °   ,  (V-A-la) 


1  -  a^ 

IT, 


2      1  -  a,  +  1  -  a. 


=   1  -  F(x^) 


-X 

=  e  °  .  (V-A-lb) 


If  we  use  these  formulas  in  Theorem  2  in  Section  III/  we  get 


and 


""2    2 
M  =   -^  X  (V-A-2) 

IT,    O 


=   —(a,  -  TT,)  (V-A-3) 

TT^     1        1 


95 


°'l  "  ""l    2 
p   =  -^ ^  x/  (V-A-4) 

TT,        O 


For  given  correlation  coefficient  p,  we  can  compute  a,  as 
a  function  of  x   from  the  formula  (V-A-4)  as 


pTT 

"l   =   — 2  "^  ^1 
o 


=   (1  +  -^)  (1  -  e   °) 


o 

IT. 


and  we  know  that  a^  =  1 (1-a,)  from  the  formulas  (V-A-la) 

and  (V-A-lb) .   From  this , 


""l 
a,   =   1  -  -^(1-  a,) 
2  ^2      ^ 


-X        X  -X    - 

O  ,     O  ,  ,         Ov  2   p 

=   e    +e(l-e    )— ^ 


2 

X 

o 


These  a,  and  a^  are  probabilities,  so  they  have  to  be  greater 
than  or  equal  to  zero  and  less  than  or  equal  to  one; 


-X 

0  £  (1  -  e    )  (1  +  -^)   £  1  (V-A-5a) 

X 

o 

—  X        X  —X 

O^e   °  +  e°(l-e  °)^  -^     <_     1  (V-A-5b) 

X 

o 

From  these  two  inequality  equations,  (V-A-5a)  and  (V-A-5b) , 
we  can  find  the  x  ranges  for  given  correlation  coefficient 
p.   To  solve  these  equations,  we  can  divide  into  two  cases. 
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one  for  positive  correlation  case  and  the  other  for  negative 
correlation  case.   If  correlation  coefficient  p  is  positive, 
then  both  equations  are  always  positive.   Thus  we  only  need 
to  find  the  x  ranges  which  makes  (V-A-5a)  and  (V-A-5b)  less 
than  1.   From  the  equation  (V-A-5a) ,  for  the  a,  case. 


-X 

a^   =   (1  -  e   °)  (1  +  -^)      1   1 

X 

o 


becomes 


/    O     ,  X  2 

p  (e   -  1)   <_  x^ 


and  let 


T  X 

^1  ==   ^o  '    y2   =   P(e   -  1) 


Because  of  the  first  derivatives  of  y,  and  y^  are  always 
positive,  we  know  that  these  two  functions  are  monotone 
increasing  functions.   Thus  we  can  find  x  ranges  which 
satisfy  y,  >^  y^  by  the  Newton  Raphson  method.   When  using 
the  Newton  Raphson  method,  let  y  =  y-i  -  y^  =  0  and  find  an 
approximate  solution,  by  approximating  exponential  series, 
which  we  can  use  as  a  starting  point.   That  is, 

2     3 

^  =  ^1-^2  =  ^o'  -  ^(^^^o^TT^^-  ^^ 


=  l^o'  -  (1-^^  ^o  ^  P   =  0 
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Then 


(1  -  -2^  -     (1  -  p  -  Yf   P  ' 
x^   =  173^ 


Starting  with  this  approximate  value  in  the  Newton  Raphson 
method,  we  can  find  x   range,  say  [x   ,x   ],  which  satisfies 
0  <_  a,  ^1.   And  for  the  a„  case. 


-x     -X      -X   ^ 

a^   =   e   °  +  e   °(l-e   O)^  ^   <   i 

^o 


becomes 


X  ^ 

p(e  °  -  1)   <   X  ^ 
—   o 


This  result  is  exactly  the  same  as  the  a,  case,  that  is,  at 

the  same  range  a,  and  a^  satisfy  constraints  a,  <_  1  and 

a^  <  1.   Thus  we  can  use  x   and  x   as  the  lower  bound  of 
2  -  ^1       ^1 

x^,x  ,  and  the  upper  bound  of  x  ,x  .   If  correlation  coeffi- 
o   £  ^^  o   u 

cient  p  is  negative,  then  the  equations  (V-A-5a)  and  (V-A-5b) 
are  always  less  than  1.   Therefore  we  need  to  consider  only 
one  constraint  which  makes  a-,  >_  0,  a2  ^  0.   From  the  a, 
equation  (V-A-5a) ,  solve  the  inequality  equation 


-X 

0  £  a^^   =   (1  -  e   °)  (1  +  -^) 

X 

o 
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-X 

since  1  -  e    is  always  positive,  we  see  that  to  satisfy 
the  inequality  1  +  — ^  should  be  positive,  i.e., 

X 

o 


>  -1 


2   - 
^o 


equivalently ,  we  have 


^o  ^ 


In  the  a-j  case,  from  equation  (V-A-5b)  , 


-X        X         -X    ^ 

0   1   a2   =   e   °  +  e  °(l-e   °)  ^  -^ 

X 


o 


or,  equivalently,  we  have 


2        /   o    T  \  2 
^o   1  -P  (e   -  1) 


As  in  the  positive  correlation  case,  we  can  find  a  starting 
point  by  approximation  to  solve  this  equation  by  Newton 
Raphson.   The  result  comes  out  as 


Xg  =  (/^  +  p)/(-  |) 


with  this  starting  point  we  find  another  bound  of  x  which 

satisfies  0  <  a^.   This  becomes  the  upper  bound  of  x  ,  x  , 
—  2  ^^  o    u 

and  from  the  a,  case,  we  have  a  constraint  x   >_  /^  which 
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becomes  the  lower  bound  of  x  ,  i.e.,  x  =  /-p.   The  lowest 

and  highest  correlations  available  for  bivariate  exponential 

pairs  in  mixture-truncation  method  are  approximately  -0.4  80  and 

0.647  respectively.  By  comparison  note  that  the  most  negative 

correlation  available  for  bivariate  exponential  pairs  with 

2 
identical  fixed  marginals  is  1  -  ^  [Moran  (1967)].   Gaver's 

(1972)  negatively  correlated  pair  has  correlations  in  the 

range  [-0.5,0].   The  table  (V-1)  shows  the  lower  and  upper 

bound  of  X   in  the  mixture-truncation  method  with  identical 
o 

marginal  exponential  and  given  correlation. 


Table  V-1:   The  lower  bound  and  upper  bound  of  x 
in  the  mixture-truncation  method  witn 
identical  exponential  marginal  distributions 

X   range  for  each  correlation 


p 

^L 

X 

u 

P 

^L 

X 

u 

0.1 

0.106 

5.832 

-0.1 

0.317 

1.984 

0.2 

0.225 

4.723 

-0.2 

0.448 

1.439 

0.3 

0.362 

3.990 

-0.3 

0.548 

1.103 

0.4 

0.527 

3.395 

-0.4 

0.633 

0.855 

0.5 

0.741 

2.842 

-0.45 

0.671 

0.751 

0.6 

1.082 

2.223 

B.   GENERATING  PROCEDURE 

The  generating  procedure  of  bivariate  exponential  random 
vector  is  almost  the  same  as  .in  the  uniform  case.   As  in 
the  uniform  case,  we  develop  three  procedures  which  are  the 
FXO  method,  the  UXO  method  and  the  TXO  method.   As  mentioned 
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in  Section  III/  all  of  these  methods  are  exactly  the  same 

except  in  how  x   is  chosen  from  the  x   range  [x  ,x  ].   The 

FXO  method  chooses  x   as  a  fixed  point  between  x   and  x  . 

This  procedure  makes  some  discontinuity  at  the  truncation 

point.   In  this  respect,  this  method  is  defective. 

The  UXO  and  TXO  methods  are  the  same  except  we  use 

different  distributions  for  x  .   In  the  UXO  method  we  use 

o 

the  uniform  distribution  and  the  triangular  distribution  for 
the  TXO  method.  These  methods  have  more  smooth  distribution 
than  the  FXO  method. 

Exponential  Mixture-Truncation  Method 

1.  (Initialization) 

i)   For  given  -0.48  <  p  <  0.64,  find  x   and  x 

2.  Define  truncation  point  x 

*  FXO  method 

^^   ^o   =  1^\    -^  ^u^ 

*  UXO  method 

i)   Generate  a  uniform  [0,1]  random  variable  U, 

ii)   x^   =   x^  +  (x^  -  x^)  *  U^ 

*  TXO  method 

i)   Generate  two  uniform  [0,1]  random  variables 

ii)   x_^   =  Xj  +  Xj_  +  X2 


where 
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X    =   (x„  +  X  )/2 

m       z  u 


^1   =   <^m-  '^t'  *  ^1 


^^2   =   '''u  -  =^m>  *  "2 


3.   Compute  parameter  values 


-X 

TT-   =   F(x^)   =   1  -  e   ° 
1        o 


"2   =   1  -  "1 


"1   =   "l*^  *  -S' 
x 
o 

a«   =   1  -  — (1  -  a,  ) 
2  TT2       1 


4.  Choose  type  for  Y 

i)   Generate  a  uniform  [0,1]  random  variable  U 

ii)   If  U  <_  7T,,  go  to  9 

5.  Y  is  an  X^ 

i)   Generate  an  exponential  random  variable  E, 

ii)   If  E,  >  X  ,  set  Y  -^  E,  and  go  to  6 
iii)   Otherwise,  return  to  5.i) 

6.  Choose  type  for  Z 

i)   Set  U  ^  ((U-  TT^)/(1-  u^)) 


lo: 


ii)   If  U  ^  1  -  Qi^,    go  to  8 


7.   Z  is  an  X2 


i)   Generate  an  exponential  random  variable  £„ 

ii)   If  E^  >  X  ,  set  Z  -«-  E«  and  go  to  11 
2    o  2      ^ 

iii)   Otherwise  return  to  7.i) 

8.  Z  is  an  X, 

i)   Generate  an  exponential  random  variable  E„ 

ii)   If  E^  £  X  ,  set  Z  ^  E2  and  go  to  11 
iii)   Otherwise  return  to  8.i) 

9.  Y  is  an  X^ 

i)   Generate  an  exponential  random  variable  E, 

ii)   If  E,  <_  X  ,  set  Y  -*-  e,  and  go  to  10 
iii)   Otherwise  return  to  9.i) 

10.  Choose  type  for  Z 
i)   Set  U  ^  U/tt 

ii)   If  U  <_  aw  go  to  8 

iii)   Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

For  the  exponential  case  it  is  possible  to  give  a  more  effi- 
cient algorithm  in  which  X,  and  X^  are  generated  exactly. 
The  algorithm  is  as  follows. 
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Efficient  Exponential  Mixture-Truncation  Method 

1.  (Initialization) 

i)   For  given  -0.48  <  p  <  0.64,  find  x   and  x 

2.  Define  the  truncation  point  x 

^      o 

*  FXO  method 

X     =    o-(X„  +  X  ) 

o     2      1  u 

*  UXO  method 

i)   Generate  a  uniform  [0,1]  random  variable  U, 

ii)   x^  =  x^  +  (^u"^e^  *  "l 

*  TXO  method 

i)   Generate  two  uniform  [0,1]  random  variables 

ii)   X   =   X,  +  x,  +  x^ 

O         K  ±  2 

where 


X   =   (x„  +  X  )/2 

m       i  u 


^1  =   (^m-^i)  *  ^1 


x^   =   (x  -  X  )  *  V^ 
z  u   m     2 


3.   Compute  parameter  values 


-X 

7Tj_  =   F(Xq)   =   1  -  e   ° 
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^2   =   ^  "  ^1 


a,   =   tTt  (1  +  -^) 


a^   =   1 -a    -    a.) 

2  TTp  1 


4.  Choose  type  for  Y 

i)   Generate  a  uniform  [0,1]  random  variable  U 

ii)   If  U  <_  TT,  ,  go  to  9 

5.  Y  is  an  X- 

i)   Generate  an  exponential  random  variable  E, 

ii)   Set  Y  ^  X  +  E, 


6.   Choose  type  for  Z 

i)   Set  U  ^  ((U-  TT^)/(1-  7T^)  ) 


ii)   If  U  <_  1  -  a^,    go  to  8 


7.   Z  is  an  X 

i)   Generate  an  exponential  random  variable  E_ 

ii)   Set  Z  -e  X  +  E-  and  go  to  11 


8 .  Z  is  an  X^ 

i)   Generate  a  uniform  [0,1]  random  variable  W^ 

ii)   Set  Z  ^  -  Ind.O  -  W2*tt-,)  and  go  to  11 

9.  Y  is  an  X^ 

i)   Generate  a  uniform  [0,1]  random  variable  W, 
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ii)       Set    Z    ^   -    Ind.O   -   W^*tt^) 


10.      Choose   type   for   Z 


i)       Set   U  ^   U/tt^ 


ii)   If  U  <_  a,  ,  go  to  8 
iii)   Otherwise,  go  to  7 


11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

Note  that  to  compute  x   and  x   in  step  1  of  both  algorithms 
we  use  subroutine  BOUND  which  is  listed  in  the  appendix  and 
the  program  with  effective  algorithm  also.   The  scatter  plot 
and  bivariate  histogram  of  resulting  random  vectors  are  shown 
in  Section  V-D. 

C.   REGRESSION  OF  Z  ON  Y  FOR  GIVEN  p 

Schmeiser  (1979)  has  used  the  regression  of  Z  on  Y  = y 

to  fix  the  parameters  in  his  bivariate  gamma  distribution. 

Consequently  we  investigate  this  for  the  mixture-truncation 

method  case.   The  regression  is  different  depending  on 

whether  Y  <  x^  or  Y  >  x_.   We  consider  two  cases  here,  one 
—  o         o 

for  fixed  x^  and  the  other  for  x  having  uniform  distribution, 
o  o 

For  fixed  x  ,  we  have 

E[ZlY  =  y,y  <_x^]   =   a^E[Xj_]  +  (l-aj_)E[X2] 

=   1  +  x^  -  a,x^(l  +  — ) 

O        1  O         IT, 
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Substituting   the  value    for    a^ , 


^O 


then  we  have 


E[Z|Y  =  y,y  ^x^]       =      1   +   x^   -    x^  (1  + -^) 

^o 


=      1   -  ^ 


^o 


And  if   y    >    X    ,    then 


E[Z|Y  =  y,y  >  x^]       =       (l-a2)E[x^]    +a2E[x2] 


1  2  ,       o 

TT,  O  TT,  2 


Substituting   the   value    for   a^, 


""l 
a         =       1    -    -—(1  -  a^  ) 
2  ^2  1 


2    x 
O 


then  we  have 


E[Z|Y  =  y,y  >  x    ]       =      1    +   ^  -^ 


TT^      X 

2      o 
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Thus  the  regression  is  constant  over  (0,x  )  and  changes  for 
y  >_  X  .   This  is  not  surprising  in  light  of  the  joint  dis- 
tribution given  in  Section  IV-D.   For  uniformly  distributed 
X  ,  the  computation  is  different  for  different  ranges  of  Y. 
If  y  <_  X  ,  then  we  have 


X 

u 
E[z|Y  =  y]   =     I   E[Z|Y  =  y,  X=x^,  ^l^o^  ^^^o^  ^^o 


If 


then  we  have 


X 

u  -X 

I   (1  -  ^)  e  °  dx 

X,  X 

I  o 


X     -X 

-X       -X.,  u  ^   o 

-  e    +  e    -  p    J    J   dx 


2  ^'^o 


""■i  ^o 


'^Z.  ~^u  ^       ,   1  ^"^u    1  ~^l 

e    -  e    +  p  ( —  e    -  —  e 

^u        ^£ 

X 

+  In  —  -  X  +  x„ ) 

x^    u    z' 


X.  <  y  <  X 

£  —  -^  —   u 
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E[Z|Y  =  y]       =  J       E[Z|Y  =  y,    X=Xq,    y>_x^]f(x^)    dx^ 


X 

u 


+      /        E[Zly  =  y,    X=x^,    y<^x^]f(x^)    dx^ 


,        -X        , 
2e"^   +    2py   -    p  (x^^  +  x^)  +  p  (^  e      ^  -  -e"^) 

u  -^ 

X 

+  p  m  ^ 
y 


If  y    >    X    ,    then  we   have 


X 

u 
E[z|Y=y]       =  I         E[z|Y  =  y,    X=x^,    y>x^]f(x^)    dx^ 


X 

U  IT-  -X^ 

1     p 


f       (1   +   -A  ^)    e      °   dx 


""l  2      o 


-X  -X 

e      ^    -    e      ^   +    p  (y   -    x^) 


By  making  x   uniformly  distributed  over  the  available  range 
of  X   for  given  correlation,  we  can  get  smoother  behavior 
for  the  regression  function. 

D.   SIMULATION  RESULTS 

We  will  show  here  the  scatter  plots  and  bivariate 
historgram  of  the  mixture-truncation  bivariate  exponential 
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random  vectors  with  correlation  p  =  0.3  and  p  =  -0.3  by 
the  FXO,  UXO  and  TXO  methods  in  Figures  (V-a) ,  (V-b)  and  (IV-c) , 
respectively.   Also  we  show  the  results  from  Gaver's  bivariate 
exponential  and  Marshall  and  Olkin's  in  Figures  (V-d)  and 
(V-e) ,  respectively. 

In  the  mixture-truncation  bivariate  exponential  distribu- 
tion, the  generated  random  vectors  by  the  FXO  method  also 
has  discontinuity  at  truncation  point  but  the  other  cases 
have  relatively  smooth  distributions. 

The  computed  correlation  printed  in  the  left  side  of  the 
bivariate  histogram  is  a  little  different  from  the  given 
correlation.   But  this  difference  can  be  assumed  as  a 
sampling  error. 

To  check  this  error  we  simulated  10  times  with  given 
correlation  p=  0.1.   From  these  simulations  we  get: 

p   =  mean  of  computed  correlation   =   0.092 

VAR[p]   =   variance  of  computed  correlation   =   0.000  8 

aCp)   =   standard  deviation  of  mean   =   0.009. 
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X 
X 


X  X 


Figure  V-al. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.   Here  x   is  fixed  at  midpoint 
between  the  lower  and  the  upper  bounds  . 
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Figure  V-a2. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
exponential  with  p  =  0.3.   Here  x   is 
fixed  at  the  midpoint  between  the  lower 
and  the  upper  bounds . 
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Figure  V-a3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.   Here  x   is  fixed  at  the 
midpoint  between  the  lower  and  the  upper 
bounds. 
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Figure   V-a4 . 


Bivariate  histogram   for    sample   of   size 
2000    from  mixture-truncation  bivariate 
exponential  with   p    =   -0.3.      Here   x      is 
fixed  at    the  midpoint  between   the  ?ower 
and  the   upper   bounds. 
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Figure  V-bl. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.   Here  x   is  taken  to  be 
uniformly  distributed  between  the  lower 
and  the  upper  bounds  . 
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Figure  V-b2 . 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
exponential  with  p  =  0.3.   Here  x   is  taken 
to  be  uniformly  distributed  between  the 
lower  and  the  upper  bounds. 
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Figure  V-b3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.   Here  x   is  taken  to  be 
uniformly  distributed  Between  the  lower 
and  the  upper  bounds  . 
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Figure  V-b4. 


Bivariate  histogram  for   sample   of   size    2000 
from  mixture-truncation  bivariate   exponential 
with    p    =   -0.3.      Here   x      is    taken   to  be 
uniformly   distributed  Between   the   lower   and 
the   upper   bounds  . 
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Figure  V-cl. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.   Here  x  has  triangular 
distribution  between  the  lower  and  the 
upper  bounds. 
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Figure  V-c2. 


Bivariate   histogram   for   sample   of    size    2000 
from  mixture-truncation  bivariate 
exponential  with    p    =   0.3.      Here   x     has 
triangular  distribution  between   tne   lower 
and  the   upper   bounds . 
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Figure  V-c3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.   Here  x  has  triangular 
distribution  between  tne  lower  and  the 
upper  bounds. 
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Figure   V-c4. 


Bivariate   historgram   for   sample  of   size 
2000    from  mixture-truncation  bi- 
variate  with   p    =   -0.3.      Here   x      has 
triangular   distribution  between   the   lower 
and   the   upper  bounds. 
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Figure  V-dl. 


Scatter  plot  for  sample  of  size  2000 
Gaver ' s  bivariate  exponential  with 
P  =  -0.3. 
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Figure   V-d2 . 


Bivariate   histogram  for   sample   of   size 
2000    from  Gaver ' s  bivariate   exponential 
with    p    =   -0.3. 
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Figure   V-el. 


Scatter  plot  for  sample  of  size  2000  from 
Marshall  and  Olkin's  bivariate  exponential 
with  p  =  0.3. 
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Figure   V-e2. 


Bivariate   histogram  for   sample  of   size   2000 
from  Marshall   and  Olkin's   bivariate 
exponential  with   p    =   0.3. 
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VI.   BIVARIATE  GAMMA  GENERATOR 

The  gamma  distribution  with  shape  parameter  r  and 
scale  parameter  X   has  the  density  function 


f  (X)   = 


r(r) 


r-1   -Ax 
X    e 


if   X  >  0 


if   X  <  0 


(VI-0-1) 


and  cumulative  distribution  function 


F(y)   =   J  j^ 


^      X^      __r-l   -AX 


(r) 


X    e    dx 


(VI-0-2a) 


In  particular,  if  the  shape  parameter  r  is  a  positive  integer 
then 


F(y) 


r-1   -Ay  , ,  . j 
1  -  I      g  ^(Ay)-* 

j=0 


(VI-0-2b) 


where  r (r)  is  Euler's  gamma  function 


r(r)   =    /   x^"-'-  e"^  dx  . 
0 


If  the  random  variable  X  has  Gamma  (r,A)  distribution  then 
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E[x)   =  f 


VAR[x]   =   -% 


when  r  =  1,  X  has  the  exponential  distribution  while  X, 
suitably  scaled,  has  an  asymptotically  normal  distribution 
as  r  ->  «■.   We  note  that  if  X  has  a  Gamma  (r,l)  distribution 
then  —  has  a  Gamma (r, A)  distribution.   So  we  may  set   =  1 

A 

as  far  as  the  generating  algorithm  is  concerned.   The  output 
from  the  generator  may  then  be  appropriately  scaled.   Two 
algorithms  for  generating  bivariate  gamma  random  vectors  are 
given  in  Section  II.   Since  they  have  complicated  computa- 
tions to  define  the  parameter  values  and  use  inverse  trans- 
formation functions,  the  mixture-truncation  method  is 
probably  preferable  to  those  methods.   To  generate  bivariate 
random  vectors  with  identical  Gamma  (r,l)  marginal  distribu- 
tions and  having  given  correlation  p  by  mixture-truncation 
method,  we  only  need  univariate  gamme  (r,l)  random  variate 
generator.   Under  the  assumption  that  the  univariate  gamma 
generator  is  available  we  developed  this  procedure.   In  fact 
we  used  Lewis  and  Robinson's  gamma  generator  as  a  univariate 
gamma  generator. 

A.   DETERMINATION  OF  PARAMETERS  IN  THE  GAMMA 
MIXTURE-TRUNCATION  METHOD 

If  the  marginal  distribution  is  gamma  with  positive 

integer  shape  parameter  r  and  scale  parameter  X  =  1,  then 
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we  know  that  F(x)  becomes  cumulative  distribution  function 
of  gamma  (r,l).   Then  the  X,  is  gamma  (r,l)  truncated  to  the 
left  of  X   and  X„  is  also  gamma  (r,l)  truncated  to  the  right 


of  X  .   Thus  we  have 
o 


X 

o 


E[Xj_]    =   y^   =     /     X  d 


F(x) 
F(x^) 


dx 


X 

1       .  °   r   -X  ^ 

TTT rTTT — r      X  e   dx 

F(x^)r(r)  ^ 


TT^rCr)^""-   ^    ^£q  (r-i)l  ""o    J^ 


X 

2   _    r  °   2  ^  F(x)       2 


VAR[X^]   =  a^^      =      J         x^  d  ^  -  .^ 


X 

o 


1     r  "  r+1   -X  ,      : 
=  V^TiFT  J        ""         e   dx  -  y^ 


=  I {(r+l)'-e   °r  y   ^^-^^^  ' X 

7Tj_r(r)^^''^-^^  •   ^    \^Q    (r+l-i)>  ""o 


-X   r+1  ,  . , . .       ^,  . 

r+l-i^^ 


2 

-  ^1 
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E[X^] 


0°  F(x)    -    F(x    ) 

2   =      I  ^  ^  — r^ 

X 

o 


F(x^) 


— =-? — r        f      X      e        dx 

2  X 

o 


-X        r 


1         r  O         V  ^i  ^~li 


TT2r(r)  ^^Q    (r-i)  i      o 


VAR[X2] 


=      a 


2         F(x)    -    F(x^) 

Turr 


I      x^    ci^-3- 


( 

7T^r(r)         ^ 


2^^^)    X 


r+1      -X    , 

X  e        dx 


,            -X      r+1         (r+1) !         r+l-i 
1        p        o      Y      ^^ —  X 

Tr2r(r)  ^^  ^£q    (r+l-i)  1      *^ 


]     -     M 


Let 


A     =      e 


-X        r 

o 


X 


r-i 


i=0 


B    =    e"''^  T    .!^r^!:.  x^^^-^ 


i=0 


(r+l-i) !      o 
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and  substitute,  then  we  have 


VARIX^l   =  .^      =   -^  -    2^    2  . 

2        ■^o 


If  we  use  these  formulae  in  Theorem  2  in  Section  III  then 
we  have 

(it  rl  -  A)  ^ 

M   =   — = 


2 

TT^Tr2(r(r))  r 


and 


a   -  TT 
^2 


Thus 


2 

(a^  -  Uj^)  (TT2r!  -  A) 

r  TT^  Tr2^(r(r))^ 


(VI-A-1) 


For  given  p,  x   and  r, 


2 

p  r!  (r-1)  !  TT,  TT^ 

a,   =   2^= —   +  TT,  (VI-A-2a) 

(TT^r!  -  A) 
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and 


TT   -  7T  (1  -  a,  ) 
a^   =  — ^  .  (VI-A-2b) 

2  TT^ 


The  value  of  (VI-A-2a)  and  (VI-A-2b)  are  determined  with 

given  p,  x   and  r.   Thus  for  given  marginal  distribution  and 

correlation  coefficient  p  we  can  find  x_  ranges  as  before. 

o     ^ 

For  simplicity,  we  will  consider  them  with  marginal  distri- 
bution gamma  (2,1).   From  the  formula  of  (IV-0-2b) , 


-X         -X 

T,     I  \  T  O  O 

^1  ^  ^^^o^   "   1  -  e    -  ^o  ® 


-X         -X 

7T^   =   1-TT,   =e    +xe 
2  1  o 


and 


A  =   e    (X  ^  +  2x   +  2)  . 
o      o 


Then  (VI-A-2a)  and  (VI-A-2b)  become 

2   P   TT    TT 

"l  =    ,  -2x *  "1  (VI-A-3a) 

o 

2   P   TT      TT 

"2  =   "2  ■"   ,   -2x    •  (VI-A-3b) 

4      o 
X   e 
o 

These  a   and  a   are  probabilities,  so  they  have  to  be  greater 
than  or  equal  to  zero  and  less  than  or  equal  to  one; 
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2   p   7T    7T 

X   e 
o 

2   p   7T^     TT 

°   i  "2   =   ^2  ^    ,   -2x    1  1 

o 


From  these  two  inequality  equations  we  can  find  the  x 

ranges  for  given  correlation  coefficient  p.   If  p  >_  0/  then 

a,  and  a„  are  always  positive.   Thus  we  only  need  to  find  the 

X  ranges  which  makes  (VI-A-3a)  and  (VI-A-3b)  less  than  one. 
o 

From  equation  (VI-A-3a) ,  the  a,  case,  and  from  equation 

(VI -A- 3b) ,  the  a^  case,  we  have  exactly  the  same  inequality 

equation; 


2pe  °(l  +  x  )   <   X  ^  +  2p  (1  +  x  )^       (VI-A-4) 
o   —   o  o 


If  p  £  0,  then  a,  and  a   are  always  less  than  one.   Thus  we 

only  need  to  find  the  x  ranges  which  makes  (VI-A-3a)  and 

o 

(VI-A-3b)  greater  than  zero.   From  the  a,  equation,  we  have 


2 


/I^d  +  X  )   <   X  (VI-A-5a) 

o   —   o 


where 


p*  =   P 


and  from  the  o.^   equation,  we  have 
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where 


/2F*  e  °   <   X   +  /27*(1  +  X  )       (VI-A-5b) 
—   o  o 


p*    = 


Note  that  if  we  let  y,  for  the  left  side  equation  of  inequality 
sign  and  ^^    for  the  right  side  equation  of  inequality  equa- 
tions (VI-A-4) ,     (VI-A-5a)  and  (VI-A-5b) ,  then  we  know  that  all 
y^  and  y^  are  monotone  increasing  functions.   Thus  as  in  the 
exponential  case  we  can  find  the  range  of  x  by  the  Newton 
Raphson  method.   For  the  positive  cirrelation  case,  equation 
(VI-A-4)  will  give  the  x„  and  x   and  for  the  negative  corre- 
lation  case,  equation  (VI-A-5a)  gives  x   and  equation 
(VI-A-5b)  gives  x  .   By  using  subroutine  GBOUND  which  is 
listed  in  the  appendix,  we  can  find  x„  and  x  .   The  tables 
(Vl-a)  and  (Vl-b)  show  the  results  from  subroutine  GBOUND. 
Unfortunately  we  have  limitations  on  the  possible  range  of 
correlations;  we  can  generate  bivariate  random  vectors  with 
correlations  in  approximately  the  range  - .  5  <_  p  <^  .  6  for  the 
Gamma  (2,1)  marginal  distribution.   In  the  same  way  as  the 
above,  the  values  x  and  x   can  be  computed  with  increasing 
complexity  for  integer  values  of  r  greater  than  2.   In 
principle  it  is  also  possible,  using  numerical  integration, 
to  compute  allowable  values  for  non-integer  values  of  r.   The 
computation  is  not  as  complicated  as  those  required  for 
Schmeiser's  bivariate  gamma.   Moreover  the  mixture-truncation 
method  does  not  require  that  the  user  be  able  to  compute 
the  inverse  gamma  distribution  function. 
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Table  Vl-a: 


X   ranges  with  Gamma  (2^1) 
marginal  distribution 
and  given  p 


p 

^L 

X 

u 

P 

^ 

X 

u 

0.1 

0.41 

7.38 

-0.1 

0.93 

3.41 

0.2 

0.65 

6.23 

-0.2 

1.18 

2.75 

0.3 

0.89 

5.42 

-0.3 

1.35 

2.34 

0.4 

1.15 

4.74 

-0.4 

1.50 

2.03 

0.5 

1.47 

4.09 

-0.5 

1.62 

1.79 

0.6 

1.95 

3.32 

Table  Vl-b:   x   ranges  with  Gamma  (3^1) 
marginal  distribution 
and  given  p 


p 

^ 

x 

u 

P 

^L 

X 

u 

0.1 

0.82 

8.99 

-0.1 

1.64 

4.73 

0.2 

1.17 

7.72 

-0.2 

1.37 

3.98 

0.3 

1.49 

6.82 

-0.3 

2.21 

3.5 

0.4 

1.83 

6.05 

-0.4 

2.4 

3.15 

0.5 

2.24 

5.31 

-0.5 

2.56 

2.89 

0.6 

2.84 

4.42 
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B.   GENERATING  PROCEDURE 

The  generating  procedure  of  bivariate  gamma  random  vec- 
tors is  almost  the  same  as  the  uniform  and  exponential  cases. 
We  also  developed  here  three  methods  which  are  the  FXO 
method,  the  UXO  method  and  the  TXO  method.   As  in  the  uniform 
and  exponential  cases,  these  three  methods  are  the  same 
except  in  the  choice  of  x   from  the  x  range  [x,,x  ].   The 
FXO  method  chooses  x   as  a  fixed  point  from  the  x   range  and 
uses  it  during  the  entire  routine  while  the  UXO  and  TXO 
methods  choose  another  x   in  every  routine.   From  experience, 
the  midpoint  of  the  x  range  gives  the  best  result  of  the 
FXO  method.   In  the  UXO  method  and  the  TXO  method,  we  assume 
that  X  has  the  uniform  distribution  and  triangular  distri- 
bution, respectively,  over  [x  ,x  ].   This  assumption  removes 
some  discontinuity  which  occurs  at  the  truncation  point  in 
the  FXO  method. 

Gamma  Mixture-Truncation  Method 
(for  integer  shape  parameters) 

1.  (Initialization) 

i)   For  given  allowable  correlation  coefficient  p , 

find  X,  and  x 
£       u 

2.  Define  truncation  point  x 
*      FXO  method 

i)    X    =   77(X„   +  X  ) 

o     2     i  u 
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*  UXO  method 

i)   Generate  a  uniform  [0,1]  random  variable  U, 
ii)   x_   =   X,  +  (x  -  X  )  *  U, 

O         I  U     £        1 

*  TXO  method 

i)   Generate  two  uniform  [0,1]  random  variables, 

ii)   X   =   X   +  X,  +  x-^ 
o       £     1     2 

where 


X   =   (x   +  X  )/2.0  , 
m       £     u      ' 


o 


''l  =   <==m-^4'*^l   ' 


^2   =   <==u-V*^2   • 


3.   Compute  parameter  values 

-X 

r-1  e    X  -' 

'^1   =   ^(-o)   =   1  -  .^   -  3 

3  =  0     -" 


where 


^2   ~   ^  "  ^1  ' 


p  rl  (r-1)  !  TT^  TT2^ 
a    =   ^ +  ^    ^ 

^  (7T2r!  -  A)^         ^ 

""l 
a^   =   1  -  — (1  -  a,  )  , 
2  TT2       1   ' 
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-X   r 


Or     r      r-1 


o   r      r 


i=0  ('^-i' 


4.   Choose  type  for  Y 

i)   Generate  a  uniform  [0/1]  random  variable  U 


ii)   If  U  <_  TT  ,  go  to  9 


5.   Y  is  an  X^ 

i)   Generate  a  gamma (r,l)  random  variable  G, 

ii)   If  G,  >  X  ,  set  Y  -*-  G,  and  go  to  6 
1    o  1     ^ 

iii)   Otherwise  return  to  5.i) 


6.  Choose  type  for  Z 

i)   Set  U  ^  ((U-  TT^)/(1  -  TT^)) 

ii)   If  U  ^  1  -  a^,    go  to  8 

7.  Z  is  an  X^ 

i)   Generate  a  gamma  (r^l)  random  variable  G^ 

ii)   If  G^  >  X  ,  set  Z  -«-  G„  and  go  to  11 
2     o  2 

iii)   Otherwise  return  to  7.i) 

8.  Z  is  an  X- 

i)   Generate  a  gamma  (r^l)  random  variable  G- 

ii)   If  G^  <  X  ,  set  Z  -*-  G^  and  go  to  11 
2  —   O  2 

iii)   Otherwise  return  to  8.i) 
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9.   Y  is  an  X^ 

i)   Generate  a  gamma  (r,l)  random  variable  G, 
ii)   If  G,  <_  X  ,  set  Y  -^  G,  and  go  to  10 
iii)   Otherwise  return  to  9.i) 

10.  Choose  type  for  Z 
i)   Set  U  ^  U/ti, 

ii)   If  U  <_  a,  ,  go  to  8 
iii)   Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

For  the  non- integer  shape  parameter  case,  step  1  and  step 
3  need  more  complicated  computation  procedures  but  here  we 
are  only  concerned  for  integer  shape  parameter  cases.   In 
step  1  we  used  the  subroutine  GBOUND  which  is  listed  in  the 
appendix.   The  scatter  plots  and  bivariate  histograms  of  the 
resulting  random  vectors  from  this  algorithm  are  shown  in 
the  next  section. 

C.   SIMULATION  RESULTS 

The  scatter  plot  and  bivariate  histograms  for  random 
vectors  of  the  mixture-truncation  bivariate  gamma  variables 
with  marginal  gamma  (2,1)  distribution  and  given  correlation 
(0.3  and  -0.3)  are  given  here.   As  in  the  uniform  and 
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exponential  cases,  the  UXO,  TXO  and  FXO  methods  are  used 

and  are  shown  in  Figures  (VJ-a) ,  (Vl-b)  and  (VI-c) ,  respectively. 

The  FXO  method  has  some  discontinuity  at  the  truncation 
point,  but  the  UXO  and  TXO  methods  appear  to  give  a  relatively 
smooth  continuous  distribution.   The  computed  correlation  in 
subroutine  BIVHST  (Bivariate  histogram)  is  a  little  differ- 
ent from  the  given  correlation.   But  we  can  assume  this 
difference  is  a  result  of  sampling  error. 
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004 


Figure   Vl-al. 


Scatter  plot  for  sample  of  size  2000 
from  mixture-truncation  bivariate  Gamma 
distribution  with  p  =  0.3.   Here  x  has 
uniform  distribution  over  (x„,x  ).  ° 
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Figure  VI-a2. 


Bivariate   histogram   for   sample  of   size    2000 
from  mixture-truncation  bivariate   Gamma 
distribution  with   p    =   0.3.      Here   x      has 
uniform  distribution  over    (x    ,x    ).  ° 
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Figure  VI-a3. 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  Gamma 


distribution  with 


p  =  -0.3. 


Here  x  has 

uniform  distribution  over  (x  ,x  ). 
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Figure   VI-a4. 


Bivariate   histogram   for   sample  of   size 

2000    from  mixture-truncation  bivariate 

Gamma   distribution  with   p    =   -0.3.      Here 

x_   has    uniform  distribution   over    (x    ,x   ) 
o  £  u 
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Figure   Vl-bl. 


Scatter  plot  for  sample  of  size  2000 
from  mixture  truncation  bivariate  Gamma 
distribution  with  p  =  0.3.   Here  x  has 
triangular  distribution  over  (x,  ,x*^). 
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Figure  VI-b2. 


Bivariate  histogram  for   sample  of   size 
2000    from  mixture-truncation  bivariate 
Gamma  distribution  with    p    =   0.3.      Here 
Xq   has   triangular   distribution   over 
(x^,x^). 
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Figure  VI-b3 


Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  Gamma 
distribution  with  p  =  -0.3.   Here  x  has 
triangular  distribution  over  (x  ,x  ;. 
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Figure  VI-b4. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma   distribution  with   p    =   -0.3.      Here 
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Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  Gamma 
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Figure  VI-c2. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma  distribution  with  p  =  0.3.  Here 
X      is    fixed  at   the  midpoint   of   x      and  x 
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Scatter  plot  for  sample  of  size  2000  from 
mixture  truncation  bivariate  Gamma 
distribution  with  p  =  -0.3.   Here  x_  is 
fixed  at  the  midpoint  of  x   and  x 
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Figure  VI-c4. 


Bivariate   histogram   for    sample  of   size 
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VII.   CONCLUSION 

The  mixture-truncation  method  is  a  general  method  which 
can  generate  bivariate  random  vectors  having  any  theoretical 
marginal  distribution  and  allowable  correlation.   The 
generating  procedure  is  very  simple  and  doesn't  need  much 
computation  for  defining  parameter  values.   In  this  respect, 
the  mixture-truncation  method  is  a  very  attractive  method  for 
generating  bivariate  random  vectors.   A  price  is  paid  for 
this  simplicity  and  generality  in  that  the  Frechet  bounds  of 
correlation  for  the  bivariate  distributions  specified  by  the 
marginal  distribution  given  by  Moran  (19  67)  are  not  always 
attained.   Also  there  is  some  discontinuity  in  the  bivariate 
distribution.   However  this  discontinuity  can  be  decreased 
by  giving  some  distribution  to  the  truncation  point  over 
its  range  for  given  p .   Thus  the  mixture-truncation  method 
is  very  attractive  for  simulation  studies  involving  only 
partly  specified  dependency  structures.   The  mixture- 
truncation  method  may  be  extended  to  generate  bivariate 
random  vectors  having  negative  values.   Another  extension 
may  be  made  to  use  grade  correlation  or  rank  correlation 
which  are  invariant  under  transformation  instead  of  using 
the  product  moment  correlation. 
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APPENDIX:   PROGRAM  LISTINGS 

r 
C 

c 

C  THE    PROGRAM    FOR    GENERATING    3IVARIATE    RAOGM    V=CTCRS(Y,Z) 

C  HAVING    IDENTICAL    MARGINAL    DI5TRI8UTICN    BY    MI  X"UR  E-"^R'JNCA 

C  C-TICN    METHOD. 

C  THE    CORRELATION    LIMTtc    ARE    AS    FOLLOWS    FOR    GIVEN    MARGINAL 

C  DISTPieU'ION. 

C  l.LNIFORM;    -0.75    <    CR    <    0.75 

C  2. EXPONENTIAL  :  -0.^8  <  CP  <  0.64 

C  i.G^MMA(2,i ) ;  -.55  <  CR  <  0.64 

C  GENERATING  PROCEDURE 

C  C^LL  BUNP(CR,NV,Y,Z,!S) 

C  CALL  BEXP(CR,NV,Yt  ZtIS) 

C  CALL  BGAM2(CR,NV, YtZtIS) 

C  WHERE 

C  CR;  GIVEN  CORRELATICN. 

C  NV;  NUMBER  OF  RANDOM  VECTORS  BE  GENERATING. 

C  Y;  the  FIRST  ELEMENT  OF  VECTO'5  WITH  DIMENSION  NV. 

C  Z;  THE  SECOND  ELEMENT  C^  VECTOR  WITH  DIMENSION  MV . 

C  IS;  INITIAL  SEED  FQR  UNIVARIATE  GENERATOR. 

C 

c 
c 
c 
c 
c 
c 


SUBROUTINE  8UNF( CR t NV t Y , Z, I J ) 


CIMENSION  Y(NV),Z{NV) 


I 


IF(CR.GT.0.75  .OR.  CR.LT.-0.75)  GC  TO  70 
1=  (CR.LT.0.0)  GO  ^0  5 
XL=C.5-( ( (9.0-12.0*CR)**0.5)/6.0) 
XU=1.0-XL 
GC  TO  10 
5  XL=(  (-1.0*CR) /3  .0**0.5 
XU=I.O-XL 
10  CO  100  I=1,NV 

Call  RANDOM(IL,Utl) 

yrj—  i(|-f(  X«J  —  XL  )  'U 

CALL    UPJOB(CR,XO,  API,  iP2,  PI  El  t  PI  E2»  BETA) 

C CHOJSE    TYPE    FOR    Y 

CALL  RANDOMi  IU,U,^  ) 
IF(L.LE.PIE1  )  GO  TO  50 
L=(U-PIE  l)/(  1.0-PIEl) 
AP=1.C-AP2 
C    Y  PRON  X2 

2  0  CALL  RANDOM ( I L,W,:  ) 
Yd  )=XO  +  W*(  l.C-XO) 
C    CHCJSE  TYPE  f=OR  Z 

IF(U.LE.AP)  GC  TO  40 
C...Z  FRON  X2 

30    CALL    RANDCM( IL,W,1) 
2(1  )=XO  +  W*(1.0-XQ) 
GO    TC    100 
C...Z    FROM    XI 

40    CALL    RANDOM( IL,W, 1) 
Z(I )=XOvW 
GO    TC    100 
C....Y    rPCM    XI 
50    U=U/FIE1 

60     CALL    RAN00M( IU,Wtl  ) 
Y(  I  )  =  XO=^W 

IF(U.LE.AP1  )     GO    TO    40 
GO    TO    30 
100    CONTINUE 
RE^LRN 
70    WRITE(6,8G)CR 
80    FGRNAT{5X, 'CHECK    CORRELATION    LIMIT, CO    NOT    ALLOW    GIVEN, 
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CCIPRELATICN 
STOP 
ENC 


SUBR 
•UNIFHR 

•FIND     I 

PIEl 
PIE2 

BETA 
RrTL' 
ENC 


GUTIN 

NITIA 
(CR  +  3 
1.0-< 
=  (  1.0 
=  (1.0 
=  AP1- 
RN 


IN    THIS    DRCGRAM    FCR    CR=»t'=7.3) 

UPJOB(CR,XO,APl,AP2,P!El,PIE2tBE"r 
VALUE 


L    VALUE 
.0*XC**2)/(3.0*XG) 
1.0-AP1)*(XQ/ (1.0-XO) ) 
-AP2)/(1.C-AP1+1.0-AP2) 
-4P1)/(1.0-AP1+1.0-AP2 ) 
(1.0-AP2) 


SIBPOUTINE  BEXP(CR»NV,Y,Z,IS) 


Cli^ENSlON    Y(NVJ  tZ(NV) 


C Ch 


C  — 


-Y 

2C 


CH 

-Z 

30 


-Z 

AO 


-Y 
50 
60 


100 

70 
80 


IX=IS 

!U=IS 

IF(CR 

CALL 

CO     10 

CALL 

Xa  =  )(L 

CALL 

casE 

CALL 

IF(L. 

L=(L- 

AP=1. 

FRGMX 

CALL 

V(  I)  = 

COSE 

IF(L. 

FRGN 

CALL 

Z(I)  = 

GC    TC 

FROV 

CALL 

2(1  )  = 

GC    TC 

FRQV 

L  =  U/P 

CALL 

Y(I)  = 

If=(L. 

GO    TC 

CGNT! 

R5TLR 

URITE 

FQRVA 

CGRPE 

STCP 

END 


*0.8 

*0.3+l 

.GT.O. 

BGUNC( 

C    1  =  1, 

RANDOM 

+(XU-X 

EPJGB( 

TYPE    F 

RANOQ^< 

LE.P^E 

DIED/ 

0-AP2 

2 

EXPGN( 

X  +  XG 

TYPE    F 

LE.AP) 

X2 

EXPON( 

X+XG 

100 
XI 

RANDOM 
-1.0+A 

100 
XI 
lEl 

RANDOM 
-1.0*A 
LE.APl 

30 
NUE 
N 

(6,8C) 
T(5Xt  • 
LATICN 


23456 

64    .OR.    CR.LT.-0.48)     GC    TO    70 

CR  tXSl »XS2tXL,XU) 

NV 

(  IU,U,1) 

L)*U 

CR  ,XG,ADl,AP2tPIEl,PI E2,EETA) 

CR    Y 

(ILItUtl) 

1)    GO    TO    50 

(1  .0-PIEl) 


IX,X»1) 


GR    Z 
GC    TO    40 

IX,X,1) 


(ILjUt 1) 
LGG(1.0-U*PIE1  ) 


(  IL,U,1) 

L0G( 1.0-U*PIE1) 

)    GC    TO    40 


CR 

CHECK    CORRELATION    LIMIT, CC    NOT    ALLCW    GIVEN, 
IN   THIS    PROGRAM    FCP    CP=',P7.3) 


SUBROUTINE    EP  JQB(  C  R  ,  XG,  API ,  AP  2  ,  P  I  El ,  PI  E2,  BE"r  A  } 
XPCNENTIAL 
PIE1=1.0-(1.0/EXP(XC)) 
FIE2=1.0-PIE1 
AP1=PI£1*(1.0  +  (CR/XC*'!=2)) 
AP2=1.0-( (PIE1/P!E2)*( 1.0-APl ) ) 
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B5TA=AP1-(1.0-AP2) 

P5TIRN 

ENC 


10 


10 


20 
25 


30 


SUBROUTINE  BOUND ( RQ  ,XS 1  ,XS2 ,X L» XU  ) 

IF{PC.LT.O.O)  GO  TO  10 

AA=(1.0-RC-((R0**2)*(5.0/12.0n  )**0.5 

XS1=( (1.0-0.5*RO)-AA)/(RO/3.0) 

XS2  =  (  (1  .0-0.5*R0)+A/!)/(R0/3.0) 

CALL    RPSON( XSltXXtRC) 

XL=  )<  X 

CALL    RPS0N(XS2,XX,RC) 

XU=XX 

PETLRN 

AR=ABS(RO) 

XL=AR**0. 5 

XS1=C.0 

XS2=((AR**0.5  )-AR)/(AR*0,5) 

CALL    RPS0N(XS2»XX,RC) 

XU  =  XX 

PETIRN 

ENC 


RPSGN(XStXX,PO) 


GO    TO 
FV,DV, 


20 

PC) 


FV,OV,RO) 


EE)     GO    TO       30 


SUBROUTINE 

EE=C.0001 

Y1  =  >S 

N  =  l 

IF(PG.LT.O.O) 

CALL     FUNP(Ylt 

GO    TO    25 

CALL    FUNNCYl, 

>2=V1-(FV/DV) 

IF( (ABS(Y1-Y2)).LE 

Y1=Y2 

GO    TO    10 

XX  =  Y2 

R=TLPN 

END 


SUBROUTINE    FUNP( Yl ♦ FV, DV, RO ) 
FV= ( Y 1**2 )-(RC*( EXP (Yl)-l.O)) 
CV=(2.0*Y1)-(R0*EXP(Y1  )  ) 

PETLRN 
END 


SUBROUTINE  FUNN( Yl , FV, DV, RO J 

RR=  ^  ES ( RO  ) 

FV  =  (Y1**2  )-PR*(  (EXP  (Yl)-1,0)**2) 

CV=2.0*Yl-( 2.C*RR*EXP(2.0*Y1)  )  +  ( 2  .0*RR*FXP ( Y 1  )) 

PETLRN 
END 


SUBROUTINE  BGAM2 ( C R tNV  ,Y , 2, IS ) 


DIMENSION  Y(NV) ,Z(NV) 

IX=IS*0.8 

IU=IS*0. 3  +  123456 

IF(CR.GT.0.64  .OR.  CR.LT.-0.55) 

CALL  GP0UND{CRtXSl,XS2,XL,XU) 

CO  100  I=1,NV 

CALL  RAN0CM(ILtU,l) 

XG=XL+( XU-XL) ^U 


GO  TO  70 
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c- 

c- 

c- 

c- 

c- 
c- 


50 
1) 


xa 

2C 


40 


CALL    GPJOB(CP  ,XG,APl,AP2,PI51tPie2tEETA) 
— CI-CaSE    TYPE    FOR    Y 

CALL  RAN0CM(IL,U,1 ) 

IP(L.LE.PIEl)  GO  TQ 

U=(L-PIE1)/(1,0-PIE 

AP=1.0-AP? 

FRCNX2 

CALL   GAMA(2.C,IX, 

IF(X.LE.XC)  GC  TO 

Y( I  )=X 
— CHCaSE    TYPE    FOR    Z 

IF(U.LE.AP)  GC  TO 

FROV  X2 

CALL  GAMA (2.0,!X,X,1 

IF(X.LE.XC)  GC  TO  30 

z(n=x 

GO  TO  100 

FRCN  XI 

CALL  GAMA(2 

IF(X.GT.XC) 

Z( l)'X 

GO  TG  100 

FRON  XI 

L=U/PIE1 

CALL  GAMA(2 

IFCX.GT.XO) 

V(I  )  =  X 

IF(L.LE.AP1} 

GG  TC  30 

CONTINUE 

RETURN 

V*Rn5(6t80)CR 

FGRMAT(5X, 'CHECK    CORRELATION    LIMIT, CO    NOT    ALLCW    GIVEN 


—  Y 

20 


— Z 

30 


—  Z 

40 


— Y 
50 
60 


100 

70 
80 


) 


0,IX,X,1) 
GC    TO    4C 


0,IX,X,1) 

GC    TO    6C 

GO    TO    40 


CCORRELATICN 
STOP 

END 


IN    THIS    PPOGRAM    FOR    CR=SF7.3) 


sue 

GANMA 
EX  = 
FIE 
FIE 
AD  = 
API 
AP2 
BET 
PET 
END 


FCUTINE  GPJ0B(CR,X0,APl,AP2,PIEl,PIE2,BE"rA  ) 

(2,1  ) 

1.0/EXP(X0) 

1=1.0-EX-XC*EX 

2=1.0-PIE1 

(X0**4)*(1.0/EXP(2.0*X0) ) 

=  PIEl+((  2,C*CR*PIE1*PIE2**2)/AC) 

=  1.0-(  (PIE1/PI»=2)*(1.0-AP1)  ) 

A=AP1-(1.0-AP2) 

LRN 


10 


SUBROUTINE  GBCUND  (  RC  ,X  S  1  ,XS  2,  XL,  XL' ) 

IF(RO.LT.O.O)  GO  TO  10 

Al=((4.0*PG)/3,0)+SQRT(  ((4,0*RO**2  )/9.0)+(4.0*RC)) 

A2-2.0*( 1.0-(P0/3,0)) 

XS1=A1/A2 

XS1=XS1+1 .0 

>S2=>S1+1C.0 

CALL  GRPSQN(XS1,XX,R0) 

>L  =  XX 

CALL  GRPSQN(XS2,XX,RC) 

>U  =  XX 

RETURN 

AR=2.0*ABS(RC) 

XL=(SGRT( ARI  +  SQRT( AR  +  SQPT(AR)*4.0)  )/2.0 

XS1=C.0 

Al=SQRT((SQRT(AR)/6.0)-KR0/9,0))-(SCRT(AR)/6.C) 

A2=SGRT(AR)/12.0 

>S2=A1/A2 

CALL  GRPS0N(XS2,XX,RQ) 

XIj  =  XX 

RETURN 
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END 
C 
C 

SU8FCUTTNE    GRPSON ( XS,XX ,R0 ) 

E5=C.001 

V1  =  )(S 
10     IF(fiC,LT.O,0)    GO    TO    20 

CALL  GFUNP(Y1  ,FV,OV,RG) 

GC    TC    25 
20    C/iLL    GFUNN(Y1  ,FV,DVtRC) 
25    >2=^1-(FV/DV) 

IF(  (ABS(Y1-Y2n  .LE.EE)    GOTO       30 

Y1  =  Y2 

GO    TO    10 
30    XX=Y2 

RETURN 

END 
C 
C 

SUERQUTINE    GFLNP ( Yl , FV, DV, RO ) 

F1=(Y1**4)+(2.0*RO*(1.0+Y1)**2.0) 

F2=2.0*R0*EXP(Y1)*(1.0+Y1) 

FV=F1-F2 

CF1  =  (4.0*Y1**3)-K4.C*RG*(  1.0+Yl)  ) 

CF2=2.0*RC*EXP(YI)*(2.0+Y1) 

CV=CF1-DF2 

FETIRN 

END 
C 
C 

SUBPCU'^INE    GFUNN(Y1  ,FV,CV,Ra) 

TR  =  SQRT(2.0*A3S(R0)  ) 

F 1=(Y 1**2.0 )+(TR*( 1.0+Yl) ) 

F?=TR*EXP(Y1) 

FV=F1-F2 

CF1=(2.0*Y1)+TR 

CF2=TR*5XP{Y1  ) 

CV=CF1-0F2 

PETLRN 

END 
C 
C 

c 
c 

c 

SUBROUTINE   BIVHST   (X,  Y,  N,  WORK) 
C 

IMPLICIT  REAL*8  (D) 

INTEGER+2   CELL,  ^^AX,  SCALE,  ONE,  CNE5 

INTEGER     RUNS,  WN,  SN,  WEIGHT,  U,  FLAG 

REAL        D,  DELTA,  OELTAX 
C 

DIMENSION  X(N),  Y(N),  WCRK(N) 

DIMENSION  XM0V(6),  YMGM(6),  OWCRK(e),  CELL(32,32) 

DIMENSION  KEY(16) ,XLABEL(8) 

DATA      NOUT/6/,  ONE/1/,  0NE5/15/ 

IF(N    ,GT.    15)    GO    TO    20 

VsRITE(NOUT,10)    N 
10    FGRWAT( •1***BIVHST***       TOO    FEW    SAMPLE    POINTS.    N    =',IS) 

PETLRN 
C 

20    AN    =    N 

AN2    =    AN    *    AN 

NV    =    N    -    1 
C CRDER    X-SAMPLE    AND    FIND    RANGES 

CALL    SORTON    (>,1,N,Y) 

XRANGE    =    X(N)    -    X(l ) 

VMA>    =    Y( 1) 

YMIN    =    Y(  1) 

CO    30    1=2, N 

IF<Y(I)     .GT.     YMAX)     YMAX    =    Yd) 
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C 


IF(Y(n  .LT.  YMIN)  YMIN  =  Yd) 
30  CONTINUE 

YR/iNGE  =  YMAX  -  YMIN 

IP(XRANGE  *  YRANGE  .QT.  0.)  GO  ^0    5C 

V.RITE(NGUT,40) 
40  FORNAK •1***BIVHST***   X  OR  Y  SAMPLE  HAS  ZERO  RANGE') 

RETURN 
C 

C 2ERC  OUT  CATA  ARRAY 

5C  DO  7C  1=1  ,32 

CO    ec   J=l  ,32 

CELL(I,J)     =    0 
60    CONTINUE 
70    CONTINUE 
C FINC    SCALE    FACTORS 

AX    =    31.999    /    XRANGE 

ex    =    1.0001    -    AX    *    X(l) 

CELTAX  =  XRANGE  /  32. 

AY    =    31.999    /    YRANGE 

PY    =    1.0001    -    AY    *    YMIN 

CELTA    =    YRANGe    /    32. 

^AX    =    0 
C ACCUMULATE    CELL    COUNTS 

DO    100    1=1, N 

IX    =    AX    *    X( I  )    +    BX 

lY    =    AY    *    Y( I  )    +    BY 

lY    =    33    -    lY 

CELLdX,     lY)     =    CELLdX,     lY)    +    ONE 

IF(CELL(IX,     lY)     .GT.    MAX)    MAX    =    CELLdX,     lY) 
100    CONTINUE 
C SCALE    FACTCR    FCR    CCLNTS 

SCALE    =    (MAX    -    ONE)     /    0NE5    +    ONE 
C FINC    SAMPLE    MCMENTS 

QWORK(l)     =    O.C 

CALL    ACCUM(X,N,nW0RK) 

CALL  VQMENT(DUORK,  XMQM) 

OWCPKd)  =  O.C 

CALL  ACCUM(Y,N,DWORK) 

C4LL  MOMENTCDWQRK.  YMOM) 
C FINC  COVARIANCE 

CSUM  =0.0 

DO  110  1=1, N 

CSUy  =  DSUM  +  (Xd)  -  XMOM(l))  *  (Yd)  -  YMHMd)) 
110  CONTINUE 

COVAR  =  DSUM  /  AN 

XDEV  =  SORT(XVOM( 2)  ) 

YOEV  =  SQRT(YN0M(2)  ) 

RHO  =  COVAR  /  (XDEV  *  YDEV) 
C KENCALL'S  TAU  COEFFICIENT 

T  =  CO 

IF{N  .GT.  500)  GO  TC  119 

CO  118  1=1, NM 

DO  114  J=I,N 

Z  =  Y(J)  -  Yd) 

IF(Z  .GT.  0.)  T  =  T  +  2.0 

IF(Z  .LT.  0.)  T  =  T  -  2.0 
114  CONTINUE 

118  CONTINUE 

T  =  T  /  (AN2  -  AN) 
C SET  UP  WORK  AREA  AND  ORDER  Y  SAMPLE 

119  CONTINUE 
VI  =  1.0 

CO  120  1=1, N 
UORK(I)  =  WI 
VI  =  WI  +  1.0 

120  CONTINUE 

CALL    SORTQN(Y,    1,     N,    WORK) 

C FINC    RANK'CRDER    S"^ATISTICS 

VN    =    0 

SN    =    C 

RUNS    =    0 

IF(Y(1)     .LT.     XI)  )     RUNS    =    1 
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FLAG  =  0 
C  =  CO 
!  =  1 
IX  =  1 
lY  =  1 
130  Z  =  ABSCFLOAT 
IF(C  .LT.  Z) 
1F(X(IX)  -  Y( 

C CURRENT  VALUE 

lAO  RUNS  =  RUNS  + 
FLAG  =  0 
lY  =  lY  ■•■  1 
1=1  +  1 
IF(IY  .GT.  N) 
GG  IC  130 

C CURRENT  VALUE 

150  RUNS  =  RUNS  + 
FLAG  =  1 
UN  =  WN  +  I 
SN  =  SN  +  WEI 
IX  =  IX  +  1 
IFCIX  .GT.  N) 
1  =  1  +  1 
GO  TO  130 

C X  AND  Y  VALUE 

160  BASE  =  X(  IX) 
N  =  1 

Jl  =  I  +  I  + 
J2  =  WEIGHTd 
K  =  I 
1=1  +  1 

C— GET  TIED  X  VA 

170  IX  =  IX  +  1 

!F(  IX  .GT.  N) 
IF(X(IX)  .NE. 
^'  =  N  +  1 
1  =  1  +  1 
Ji  =  Jl  +  1 
J2  -  J2  +  WEI 
GO  TC  170 

C FINC  TIED  Y  V 

180  lY  =  lY  +  1 

IF(IY  .GT.  N) 
IF(Y(IY)  .NE. 

1  =  1  +  1 
Jl  =s  Jl  +  I 
J2  =  J2  +  WEI 
GO  TO  180 

C WEIGHTED  AVER 

190  L  =  I  -  K  +  1 

2  =  FLO  AT (M  * 
V»N  =  WN  +  Z 

Z  =  FLQAT(M  * 
SN  =  SN  +  Z 

C AD  FCC  CORREC 

PUNS  =  RUNS  + 
1  =  1+1 
iFdX  .GT.  N) 
IF(IY  .LE.  N) 

C Y  CeSERVATICN 

200  RUNS  =  RUNS  + 
N2  =  N  +  N 
DO  210  1=1, N2 
WN  =  WN  +  I 
SN  =  SN  +  WEI 
210  CONTINUE 

C FINC  MANN-WHI 

22C  L  =  WN  -  (N  * 

C SPEARMAN'S    RA 

CSUN    =    0. 

WI    =    1.0 

DO    230    1=1, N 


(IX-TY))    /    AN 

C    =    Z 

lY))     150,     160,    140 

IN    PGCLED    SAMPLE    IS    A    Y 

FLAG 


GO    TO    ?00 

IN  PGCLED  SAMPLE  IS  AN  X 
1  -  FLAG 


GHTd  ,N) 
GO  TO  2?0 

S    ARE    TIEC.    USE    MIDRANK    CORREC^ICN 

1 

,N)     +    WEIGHT(I+l ,N} 


LUES 

GO    TO    180 
BASE)    GC    TO    180 


GHTd  ,N) 

ALUES 

QQ    TO    190 
BASE)    GO    TC    190 

GHT( I ,N) 

AGE  OF  TIED  VALUES 

Jl)  /  FLCAT(L)  +  0.5 

J2)  /  FLaAT(L)  +  0.5 

TICN    FCR    RUNS 
FL0AT(2    *    M    *    (L-MJ  ) 

GO    TO    22C 
GO    "^0    130 
S    EXHAUSTED 
1 

GHT( I,N) 


/    =LnAT(L) 


TNEY    STATISTIC 

(N    +    1))     /    2 
NK   CORRELATION    COEFFICIENT 
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CSU^^    =    DSUM    +    (WI    -    WCRKd))     **    2 

WI     =    WI    +    1.0 
230    CONTINUE 

R    =    1.    -    6,    *    OSUM    /    (AN    *    (AN2    -    1.)) 
C FINC    SAMPLE    MEDIANS 

IF(N0D(N,2)     .EQ.    0)     GO    TO    240 

N    =    1    +    N/2 

XMEC    =    X(M) 

>(MEC    =    Y(  N) 

GO    TQ    250 
240    N    =    N/2 

XMEC    =    0.5    *    {X(M)     +    X(M+1) ) 

YMEC    =    0.5    *    (Y(M)     +    Y(M+1) ) 

C FIND    NORMALIZED    STATISTICS 

250    PNCPM    =    (RUNS    -   N    -    1)     *    SQRT((AN    +    AN    -    I.)     /    (AN2    - 
CAN)  ) 

FAC    =    SQRT(12.    /     ( AN2    *    (AN    +    AN    +    1.))) 

UNORM    =    (U    -    C.5    *    AN2)    *    FAC 

VNOPM    =    (WN    -    AN2    -    0.5    *    AN)     *    FAC 

SNCRM    =    ( SN    -    AN2    -    0.5    *    AN)     *    FAC 
C REORDER    Y    SAMPLE 

CALL    SORTGN     (WORK,     1,    N,    Y) 
C 

C  V«RI7E    OUTPUT 

C 
C hEACING 

WRITE(NOUT,    3C0)    N 
300    FGRMAT( 'l*  ,20Xt'BIVARIATE    SAMPLE    SUNMARY ', lOX  ,» SAMPLE 
C    SIZE    =S  19//) 

V»RITE(N0UT,31C) 

310    FQR^'AT(13X,  •+  — •  ,7(«*»»7(»-'))t  •* +•  ) 

C 

C eCD\    OF    PLOT 

CALL    0UTPU^(C=LL(1 t    l)t    SCALE,    1,     C.) 

CALL    REPEAT 

CALL    OUTPLT(CELL( 1  ,    2),    SCALE,       2,     VMAX    -  DELTA) 

VRIT£(N0UT,32C) 
320    FORVAT(«+« ,93X, 'MEASURES    OF    ASSOCIATION') 

CALL    REPEAT 

CALL    OUTPUT(CELL( 1,    3),     SCALE,    3,     C.) 

V»RITE(NCUT,330)    COVAR 
33  0    FGRMAT( •  +  ' , 83X ,' COV ARI ANCE ' ,2 5X, IPE  14 .6 ) 

CALL    REPEAT 

V.RITE(N0UT,34C)  RHC 
340  FaRNAT( •+',83X, 'CORRELATION  COEF F IC I 5NT' , 13XF9 . 6 ) 

CALL  CUTPUT(CELL(1  ,  4),  SCALE,  4,  0.) 

URnE(N0UT,35C)  R 
350  FQPMAT(»  +  «  ,83)<, 'SPEARMAN' 'S  RANK  CCRRELA^ICN  CCEF. 
C,F9.6) 

CALL  REPEAT 

IF  (  N  .LE.  500  )  WRITE(N0UT,360)  T 
360  F0RMAT( '-f*  ,83X,'K?NDALL"S  TAU  CCEFF  ICI  ENT  '  ,  1 IX  F9  ,6  ) 

CALL  0UTPUT(CELL(1  ,  5),  SCALE,  5,  C.) 

CALL  REPEAT 

CALL  0UTPUT(CELL(1,  6),  SCALE,   6,  YVAX  -  5 .  =»  DELTA) 

V»RITE(N0UT,37C) 
370  FGRNAT(»  +  '  ,95X, 'TESTS  FOR  EQU  ID  ISTR  IBUT  ION  '  ) 

CALL  REPEAT 

CALL    CUTPUT(CELL(1  ,    7),    SCALE,    7,     0.) 

WRITE(N0UT,38C) 
380    FCPNAT('+',111X,'TEST' , 7X, • NORMAL  I  ZED' ) 

CALL    REPEAT 

WRITE(NOUT,39C) 
390    FGRVAT(«  +  '  ,109X,'  STATISTIC  ,5 X,  'STATISTIC  ) 

CALL    0UTPUT(CELL(1  ,    8),     SCALE,    8,    0.) 

V<RITE(NOUT,40C)    0 
400    FaRMAT('+' ,83X,'K0LVCGCF0V-SMIRN0V    TEST    •,F10.e) 

CALL    REPEAT 

WR  ITE(N0UT,41C)  RUNS,  RNORM 
410  FaRMAT( •+',83X,'WALD-WCLF0WITZ  RUNS  TEST'  ,  1 10  ,  5XF9 . 5  ) 

CALL  0UTPUT(CELL(1  ,  9),  SCALE,  9,  0.) 

WRITE(N0UT,42C)  U,  UNQRM 
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420  FaRMAT('  +  '  ,83>,•MAN^-WHITNEY  TEST '  ,  7X I10» 5XF9  .5 ) 

CALL  REPEAT 

V»RITE(N0UT,43CI  WN,  WNCRM 
430  F0RVAT<«+«,83X,»WlLCaXS0M  TEST' ,  lOXI  10  ,5X1=9  .  5  ) 

CALL  CUTPUT(CELL(1  tlO)  ,  SCALE,  10,  VHAX  -  9 .  =>  DELTA) 

URITE(N0UT,44C)  SN,  SNCPM 
440  FaR^'AT(  •  +  «,83>, 'STEGEL-TUKEY  TEST  •  ,  7X1 10,  5XF9  .5  ) 

CALL  REPEAT 

CALL    0UTPIJT(CELL(1  ,11),    SCALE, 11,    0.) 

CALL    REPEAT 

V<RITE(N0UT,45C) 
450    FORMAT( •  +  •  ,92>, 'UNIVARIATE    STATISTICS') 

CfiLL    OUTPUT(CELL( 1,12),    SCALE, 12,     C.) 

CALL    REPEAT 

WRITE(N0UT,46C) 
460  FCPVAT( •+«,g5X,'X  S AMPL E' , 9X, • Y  SAMPLE') 

CALL  0UTPUT(CELL(1  ,13),  SCALE, 13,  0.) 

WRITE(N0UT,47C)  XMOM(l),  YMOMCl) 
470  FORMAT ( '+ '  , 83X , ' ME  AN • , 5X1PE14.6 , 3XE 14 .6 ) 

CALL  REPEAT 

WRITE(N0UT,48C)     XMEC,    YMED 
480    FaPVAT('  +  '  ,83X,'MEDIAN»,3X1PE14.6,3XE14,6) 

CALL    0UTPUT(CELL(1,14) ,    SCALE,    14,    YNAX    -13.     *    DELTA) 

CALL    REPEAT 

WRITE(M0UT,49C)  X^^0N(2),  YMaM{2) 
490  FQRWATC '+' , 83 X, ' VARI ANCE » , IX IPE 14. 6 , 3XE14. 6 ) 

CaLL  0UTPUT(CELL(1 ,15),  SCALE, 15,  C.) 

V»RITE(NOUT,50C)    XOEV,    YDEV 
500    FGRV^TC »+• ,83X, 'STO    DE V ' , 2X1PE14. 6 , 3 XE14. 6 ) 

CALL    REPEAT 

l«RITE(N0UT,51C)  XRANGE,  YRANGE 
510  FCRVAT(»+',83X,'RANGE',4X1PE14.6,3XE14.6) 

CALL  OUTPUT(CELL( 1 ,16),  SCALE, 16,  0.) 

CALL  REPEAT 

WRITE(NOUT,52C)    XMCM(4),     YM0M(4) 
520    FaPVAT(»+«  ,83X,'SKEWNESS' ,1 X1PE14 .6 ,3XE14.6 ) 

CALL  0UTPUT(CELL(1 ,17) ,  SCALE, 17,  0.) 

V^RITE(NOUT,53C)    XMCN(6),     YMGM(6) 
5  30     FaRNAT('+',a3>,'KURTGSIS',lXlPEl4.6,3XE14.6) 

CALL    REPEAT 

CALL    OUTPLTCCELLd  ,18)  ,    SCALE,    18,     V^'AX    -17.     *    CELTA  ) 

V<RITE(N0UT,54C)     X(N),    YMAX 
540    FCRNATC •+• , 83X , ' MAX IMUM ' , 2X1PE14. 6  ,  3 XE14. 6  ) 

CALL    REPEAT 

^RITE(NOUT,55C)    X(l),    YMIN 
55C    FQRNA^('+',83X, » MI N IMUM ' , 2X IP E14.6, 3X= 14.6 ) 

Z    =    YMAX    -    17.    *    DELTA 

DC    560    J    =    19,32 

Z    =    Z    -    DELTA 

CALL    OUTPUT(CELL( 1,J),     SCALE,    J,     Z) 

CALL  REPEAT 
560  CONTINUE 

WPITE(N0UT,31C) 

XLAEELd)  =  X(l)  +  OELTAX 

FGUPD  =  4.  *  CELTAX 

CO    570    1=2,8 

XLAEELd)    =    XLABEL(I-l)    +    FCURD 
570    CONTINUE 

WRITE(N0UT,58C)(XLAeEL( I),I=l,7,2),(XLAeFL(J),J=2,8,2) 
580    F0RMATd6X,8(  '1  '  ,7X)/    12X  ,4dPE10  .3  ,  •        I  •)/    20X,4(E 

C1C.3,6X)  ) 

KEYd)    =    0 

CO    590    1=2,16 

KEYd  )    =    SCALE    *    ^    -    1 
590    CONTINUE 

V,RITE(NOUT,60C)     (  K  EY  ( I  )  ,1  =  1 , 1  6  ) 
600    FGRVAT(//50X, 'KEY' ,//•     SYMBOL    PR  INTED* , 12X , • - 

1-  =  r  =  =  T>TA 

2    V  H  I-  H»/       •+•  ,38X,  '  .  .  !  * 

34X<XVE  EEE' 

4/       •+',56X,'.  .  -  -  T  = 

5/  T  X'/    •+'  ,110X, 'T'// 

162 


6   •  NO.  OBSERVATIONS', 15, 1516) 

RETURN 

END 

FUNCTION  WEIGHT(I»  N) 

INTEGER  WEIGHT 
C 

C      THIS  FUNCTION  SUBPRCGRAM  COMPUTES  ThE  THE  VvETGHT  FUNCT 
C      ICN  FOR  THE  I-TH  RANK  IN  A  SIEGEL-TUKEY  TEST  VnITH  A  SA 
C      -MPLE  OF  N. 
C 

K  =  fODd  ,2) 

UETCHT  =  I  +  I  -  K 

IF(I  .GT.  N)  WEIGHT  =  4  *  N  -  WEIGHT  +  2  *  (1  -  K) 

RETLRN 

END 

SLBRCUTINE  OUTPUT  (CCUNTS,  SCALE,  L,  Y) 
C 

C      THIS  SUBROUTINE  OUTPUTS  A  SINGLE  LINE  OF  THE  PIVARIATE 
C      HISTOGRAM  CN  THE  LINE  PRINTER 
C      ON  THE  LINE  PRINTER 
C 

INT5GER*2   COLNTS,  SCALE,  ZERO,  ONE 

LOGICAL*!   LINE,  CHAR 

DIMENSION   LINF(64,A),  CHAR(4,  16),  C0UNTS(32),  N(16), 
CCHA(16) 

EQUIVALENCE  (CHAR(1,1),  CHA(l)) 

DATA        ONE/1/,  ZERO/0/,  NaUT/6/,  N/3*l, 3*2 , 2*3, 2, 6 
C*3,^/ 

CAT4        CHA  /•     «,•-    ',»=    •,'-.   ','=.   •,'=! 
1   •,»=*.  •,  •=+.  •♦•TX   •,•><-  N'TX- 

2S«AVT  •,'HE=  S'HE/  »,  •  HET  »,«HEXT»/ 

C 

LINES  -  1 

K  =  1 

CO  20  J=l  ,32 

INDEX  =  CCUNTS(J)  /  SCALE  +  ONE 

IFdNDEX  .EQ.  1  .AND.  CCUNTS(J)  .GT.  ZERO)  INCEX  =  2 

IP(N(INDEX)  .GT.  LINES)  LINES  =  N(INCEX) 

DO  10  1=1,4 

LINE(K,I)  =  C»-AR(  T,  INDEX) 

LIN£(K+1,I)  =  CHAR(  T, INDEX) 
10  CONTINUE 

K  =  K  +  2 
20  CONTINUE 

IF(V0D(L-2,4)  .EQ.  0)  GO  TO  40 

ENTRY  REPEAT 

WRITE{N0UT,30)  ((LINE(J,I),  J=l,64),  1=1, LINES) 
30  FCRNflT(i3x,« 1  •,64A1,«|  •/(•  +  ', 13X , 64A1) ) 

GO  TO  60 
40  V«RITE(NGUT,50)  Y,  ((LINE(J,I),  J  =  l,64),  1  =  1, LINES) 
50  FOPVAT(1X1PEI0.3,  •   *',64A1,'*'/  (  •  +  • , 13X , 64A  1 )  ) 
60  RETLRN 

END 
C 
C 

SUBROUTINE  MOMENT  (WORK,  XMOM) 
C 

C      PLRFCSEi 

C         FINO'tHE  MOMENTS  OF  THE  DATA  ARRAY  X  WHEN  NCT  ALL 
C      CF  THE  ARRAY  IS  AVAILABLE  AT  ANY  GIVEN  TIME. 
C 

C         CATA  DOINTS  ARE  PASSED  TO  THE  SIEROU^INE  "N"  AT  A  T 

C  -IMC  WITH  THE  RUNNING  TQTALS  KEPT  IN  "^HE  DTUBLE  PRECIS 

C  -ICN  WORK  AREA  VECTOR  'WORK'.   PRIOR  TO  THE  FIRST  CALL 

C  ,W0RK(1)  SHOULD  BE  SET  TO  0.0.  DATA  IS  THEN  P/5SSE0  BY 

C  THE  CALL 

C  CALL  ACCUM(X,  M,  WORK) 

C  AN  INITIAL  MEAN  ESTIVATP  IS  KEPT  IN  UnRK(4)  AFTER  THE 

C  FIRST  7  DATA  FCTNTS  HAVE  BEEN  PASSED.   HIGHER  CENTRAL 

C  NCMEN'^S  ARE  FCUND  USING  THIS  ESTIMATE.  ONCE  AT  LEAST  7 

C  X  VALUES  HAVE  BEEN  ACCUMULATED,  RESULTS  MAY  RE  CBTAINE 
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C  -0     FROM 

C  CALL    MCMENT(UCRK,    XMQM) 

C  PESLLTS    ARE    RETURNED    AS    FOLLOWS: 

C  >NGM(1)       THE    SAMPLE    MEAN 

C  >M0M(2)       THE    SAMPLE    VARIANCE    (UNEIASEO    ESTIVATOR) 

C  )(^'0M(3)       THE    SAMPLE    THIRD    CENTRAL    MOMENT     (LNEIASED 

C  ESTIMATOR) 

C         >MCM(4-)   THE  SAMPLE  COEFFICIENT  CF  SKEWNESS 

C         XN0M(5)   THE  SAMPLE  =OURTH  CENTRAL  MQMEN^  (UNBIASED 

C  ESTIMATOR) 

C         >MCM(6)   THE  SAMPLE  COEFFICIENT  CF  KURTQSIS 

C 

C 

IMPLICIT      REAL*8     (D) 

R'~AL^8  UCRK 

DIMENSION   W0RK(8),  >MCM(6}t  X(M) 
C 

TF(V»0RK{1)  .LT.  7,0D0)  RETURN 
C FIND  CORRECTION  FACTORS 

CMEAN  =  <WQRK(2)  +  WCRK(3))  /  WORK(l) 

DC  =  DMEAN  -  W0RK(4) 

CC2  =  DC  =»  DC  ♦  WORK(I) 
C DETERMINE  CORRECT  CENTRAL  MOMENTS 

DM2  =  W0RK(5)  -  002 

DM3  =  VI0RK(6)  +  WCRK(7)  -  DC  *  (3. CO  "*  W0RK(5)  -  002  - 
CCC2  ) 

DM4  =  W0RK(8)  -  DC  *  (4.00  *  (W0RK(6)  +  W0RK(7))  -  PC 
1  *(6.C0  *  WCRK(5)  -  3. DO  *  CO?)) 

C CORRECT  ESTIMATORS  IN  WORK  ARRAY 

V*0RK(4)  =  DMEAN 

^aRK(5)  =  DM2 

K0RK(6)  =  0. 

WCRK{7)  =  0. 

IF(CM3  .GT,  0.)  W0RK(6)  =  DM3 

IF(CM3  .LT.  0.)  W0RK(7)  =  DM3 

V.QRK(8)  =  DM4 
C PETLRN  STATISTICS 

XMCNd)  =  DMEAN 

CFAC  =  WORKCl  )  -  1 .DO 

XMCV(2)  =  0M2  /  DFAC 

CFAC  =  DFAC  *  (WORKCl)  -  2 . DO ) 

>MCN(3)  =  WORK(l)  *  DM3  /  DFAC 

XM0M(4)  =  XMCV(3)  /  (XMQM(2)  *  SQRT  (  XMOM( 2 )  )  ) 

CFAC  =  DFAC  *  (WORK(l)  -  3.00) 

XM0V(5)  -    (((WORK(l)  -  2. DO)  *  WQRK(l)  +  3.^0)  *  0M4 
1  -  (6.00  *  WORK(l)  -  "5. CO  *  DM2  ^    DM2/ 

2V»aRH(l))  /DFAC 

XMQV(6)  =  XMOfCS)  /  (XMGM(2)  *  XM0M(2))  -  3.0 

PETLPN 
C 
C 

c 

IF(VsQRK(l)  .GT.  6. DO)  GO  TO  210 

C ACCUMULATE  FIRST  7  DATA  POINTS 

N  =  V^ORK(  1)  +  2.00 
V»aRK(l)  =  WORK(l)  +  M 
NO  =  ^  +  V  -  1 

IF(NP  .GT.  8)  NP  =  8 
J  =  1 

CO  201  I=N,  NF 
V«ORK(I)  =  X(J) 
J  =  J  +  1 
2C1  CONTINUE 

IF(NP  .LT.  8)  RETURN 

C CETMN    INITIAL    MEAN    ESTIMATE 

CSUVP  =  0. 
CSUMM  =  0. 
CO  202  I  =  2,8 

!=(V»ORK(!)  .LT.  0.)  D5UMM  =  DSUMM  +  WGRK(!) 
IF(V»aRK(I)  .GT.  0.)  DSUMP  =  OSUMP  +  WORK(I) 
202  CONTINUE 


ENTFY  ACCUMCX,  M,  WCRK) 
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2C3 

204 


IF(J 

CO  2 
IF(X 
IF(X 
CQNl 
CMEfl 


205 


;— SET 

D2  = 
C3M 
C3P 
C4  = 
CO  2 
CIFF 
CF2 
C2  = 

I5=(C 

!F(C 

C4  == 

CONT 

IF(  J 

CC  2 

CIFF 

CP2 

C2  = 

IF(C 

IF(C 

C4  = 

CCNT 

WORK 

VQRK 

V«OPt< 

KORK 

WCRK 

UORK 

WORK 

RETU 

ACCl 

210  WCRK 
CO  2 
IFO 
IFO 
CIFF 
CF2 
V^QRK 
IF(C 
1F(C 
V^ORK 
CQMT 
PETL 
END 


206 
207 


•GT.  M) 
03  I=J,y 
(I)  .LT. 
(I)  .GT. 
INUE 

N  =  (CSUM 
UP  WORK  A 

C, 
=  0. 
'    0. 

0. 

05  1=2,8 
-    WORK (I 

'    DIFF  ♦ 
02  +  DF2 

IFF  .LT, 

IFF  .GT. 
04  +  D«=2 

INUE 
•GT.  M) 

06  I=J,M 
=  X( I)  - 

=  DIFF  * 

02  +  DF2 
IFF  .LT. 
IFF  .GT. 
+  DF2 


GO  TO  204 

C.)  nSUMM  =  DSUMM  +  X(I) 
C.)  OSUMP  =  DSUMP  +  X(I) 


N  + 
REA 


DSLMP)  /  WQRK(l) 


OSU 

CSU 

DME 

02 

D3P 

03^ 

D4 


C- 


211 


D4 
INUE 
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
(8) 
RN 

^'LLATE  DA 
(1)  =  WOP 
11  1  =  1, ^^ 
(I)  .GT. 
(I)  .LT. 

=  X{  I)  - 
=  DIPP  * 
(5)  =  WOR 
IFF  .GT. 
I^F  .LT. 
(8)  =  WOP 
INUE 
RN 


)  -  OMEAN 
CIFF 

C.)  D3N  = 
C.)  03P  = 

*  DF2 

GO  TO  207 

CMEAN 
CIFF 

0.)  03V  = 

c.)  n3P  = 

*  0'=2 
NP 

AN 


C3M 
C3P 


DIFF 
DIFF 


CF2 
DF2 


C3M 
C3P 


DIFF 
DIFF 


CF2 

CF2 


TA 
Kd) 


+  M 


C.)  W0RK(2)  =  W0RK(2)  ♦  X(I) 
0.)  W0RK(3)  =  W0RK(3)  +  X(!) 

WGRK(4) 
CIFF 

K(5)  +  OF? 

C.)  WCPK(6)  =  W0RK(6)  ■»  CIFF 
C.)  WCRK(7)  =  wnRK(7)  -*  DIFF 
K(8)  +  CF2  *  0^2 


DF2 
DF2 


SUBPGUTINE   SCRTON  (ON,  II,  JJ,  WITF) 

THIS  SUBROUTINE  SORTS  ThE  VECTOR  ON  INTO 
ASCENDING  ORDER  FRON  CN(II)  TO  ON(JJI,  CARRYING 
ALONG  THE  CORRESPONC ING  ELEf^ENTS  OF  THE  VECTOR 
WIT1-. 


INTEGER  HI  ,  HISTK, 
OIMENSION   QN(JJ), 


FIKEEP,  STKPT,  WITH 

WITH(JJ),  HISTK(16),  L0STK(16) 


10 
20 


STKPT  =  1 

LOKEEO  =  II 

hIKEEP  =  JJ 

IF  (LOKEEP  .GE. 

LO  =  LOKEEP  -  1 

FI  =  HIKEEP  +  1 

MIDCLE  =  (LOKEEP 

IL  =  LOKEEP 

IN  =  HIKEEP 

IF  (  ON(MIDOLE)  .GT 


HIKEEP)   GO  TO  100 


+  HIKEEP)  /  2 


ON(HIKEEP)  )  IF  =  MIDDLE 
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30 


50 
70 


90 


100 


110 


120 


130 


140 


IF 
IF 

IF 
IS 
T  = 
GO 

IF 

T  = 

GN 

CN 

IT 

UI 

WI 

HI 

IF 

LC 

IF 

IF 

IF 

LQ 

l-I 

ST 

LQ 

QQ 

LQ 

HI 

ST 

HI 

GO 

ST 

IF 

LC 

HI 

IF 

IF 

LQ 

LC 

IF 

IF 

T5 

LQ 

CN 

LC 

IP 

CN 

IW 

N 

IT 

CO 

V«I 

IW 

CO 

WI 

GO 

EN 


(  QN(MIDOLE) 
(  nM(HIKeEP) 
(  aN(L0K5EP) 
UE  =  LGKEEP  + 
SI  =  ONdSUE) 
■TO  50 


.LT.  ON(LQKEEP)  )  IL  =  MIDDLE 

.LT.  ON(IL)       )  IL  =  HIKEEP 

.GT.  ON(IH)      )  IH  =  LQKEFF 

MIDDLE  +  HIKEEP  -  IH  -  IL 


( 

MP 

(H 
(L 

EV 
TH 
TH 

'( 

~{ 
( 

(H 

SI 
S7 
KP 

KE 

T 

ST 

ST 

KP 

KE 

T 

KP 

(5 

KE 

KE 

( 

( 

KE 

KE 

( 

< 

MP 

(L 

( 

(L 


EV 

1 

TH 

Nl 

TH 

T 

D 


LO  •£ 

I)  = 

0)  = 
P 

(HI) 
(LO) 

HI  - 

ON(H 

LO  + 

CN(L 

LO  . 

I  -  L 

K(STK 

K(STK 

T  =  S 

EP  = 

G  110 

K(STK 

K(STK 

T  =  S 

EP  = 

C  110 

T  =  S 

TKPT 

EP  = 
EP  = 
HIKE 
LOKE 
EP  = 
EP  = 
LOKE 
aN(L 
=  CN 
LCKE 
0+1) 
LC  - 
TEMP 
0  +  1) 

LCKE 
LCKEE 
P  =  W 
40  IM 
(IW) 

IW  - 
INUE 
(IW) 
0  120 


CN 
CN 
TE 


1 

I) 

1 

0) 
LE 
CK 
PT 
PT 
TK 
LQ 


HI)  GO  TO  50 
(HI) 
(LC) 

MP 

WITH(HI ) 
WITH(LQ) 
ITEMP 

.GT.  TEST)  GO  TO  50 

.LT.  TEST)  GO  TO  70 
.  HI  )  GC  TO  30 

EEP  .LE.  HIKEEP  -  LO)  GO  TO  90 
)  ^    LOKEEP 
)  =  HI 
PT  +  1 


PT)  =  LO 
PT)  =  HIKEEP 
TKPT  +  1 
HI 


TKpT 

.EQ. 

LOST 

HIST 

EP  - 

EP  . 

LCKE 

LCKE 

EP  . 

CKEE 

(LOK 

EP 

=  CN 

1 

.LT 
=  TE 
EP  + 
P  - 
ITH( 
CVE 
=  WI 

1 


TO 


11) 

10 


GO  TC  20 


-  1 

0)  RETURN 
K(STKPT) 
K(STKPT) 

LOKEEP  .GE 
EQ .  II)  GO 
EP  -  1 
EP  +  I 

EO.  HIKEEP  )  GO  TO  ICC 
P)  .LE.  CN(LaKE£P+l)  ) 
EEP  +  1) 

(LO) 

.  CN(LC) )  GO  TO  130 

^'P 
1 

LO 
!W) 

=  It  N 
IH(IW-I) 


GO  TO  120 


=  ITEMP 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


FCP    GfNER/STING    EIVARI^TE    RANDOM    V  ECTORS  ( Y,  Z  )  ,C  ALL 
FCLLCWING    SUeRGLTIN5. 

!•    INIFORM(0,1)    MARGINAL    DI STRI BUT ICN. 

CALL    LEV^UFh'(CR,NV,Y,ZtIS);     POSITIVE    CORREL/!TION    BY 

LAWERENCE    t^HD    LEV^IS    METHOD. 

C/!LL    GAVUF^(CR,NV,Y,Z,  IS  )  ;    NEGATIVE    CORRELATION    BV 

TRANSFCRNA7I0N    CF    GAVER'S    EXPa^E^TIAL. 
2.     EXPCNENTIAL(MEAN=1.0)    MARGINAL    C  I  S  TR  I  BL'TIC^  . 

CALL    MAREXP(CR,NV,Y,Z,IS)  ;     POSITIVE    CORREL/STION    BY 

MARSHAL    ANC    OLKINS    NETHOD. 

CALL    GAVEXF(CR»NV,Y,Z, IS);    NEGATIVE    CORRELATION    BY 

GAVER'S    NETHQO. 
WHERE 

CR;    GIVEN    CCRRELATICN. 

NV:     NLMBEP    OF    RANOCr'    VECTORS    BE    GENERATING. 

Y;    THE    PIRST    ELEMENT    OF    VECTOR    WITH    CIMENSICN    NV. 

2;    THE    SECONC    ELEMEN"^    CF    VECTOR    WITH    DIMENSION    NV. 

IS;     INITIAL    RANDOM    SEED. 


10 
50 


6C 
100 

2C 
30 


SUBPCUTINE  LEWUFM{CP,NV,Y,Zt lU) 


01 
IF 
A  = 
BE 
AL 
F  = 
R  = 
CO 
CA 
Y( 
IP 
U  = 
IF 
2( 
GC 
2( 
GO 
V  = 
IF 
Z( 
GC 
Z( 
C3 
PE 

FC 

EN 


MEN 
(C 
SCR 
TA  = 
P  t- 
(1. 
1.0 
IC 
LL 
I  )  = 
(L( 
(L( 
(U. 
I  )  = 
TC 
I  »  = 

TC 

L(2 
(V. 
I)  = 

TC 
I  )  = 

NT  I 
TLR 
HE 
RNft 
TLR 
C 


SION    Y(NV) ,Z(NV),U(3) 

R.LT  .0.0)    GO    TC    20 

T(36.Q+(  12.0*CR)  +  (16.0*(CR**2.C)  ))--6.0 

A/(2.0*CP) 

2.0-(1.0/BETA) 

C-BETA)/(1.0-BETA+(ALPA*BETA) ) 

-ALPA 

0    1=1, NV 

RANDOM( IL,U,2) 

U(l) 

2)  .LE.R)    GO    TC    50 

2)-R)/( l.O-R ) 

LE.P)    GC    TO    10 

(Y(I)**BETA)*((W-P)/(1.0-P))«*(f5*BE'A) 

100 
(Y(I)**B5TA)*(U/P) 

100 
)/R 

LE.P)    GC    TO    60 
((W-P)/(  1.0-P)  )**(R*BETA) 

100 
V/P 
NUE 
N 


(6,30)CR 
T(5X,«THIS 

N 


METHOD    NOT    ALLOW    CR=»»F7.3) 


SUBROUTINE    GA VUFM( C P , NV tY, Z »I  S) 

DIMENSION    Y(NV),Z(NV),WK(1),G(100) 

COUELE    PRECISION    DXS 

CXS=IS*0.8 

IL=IS*0. 5+123^5 

IF     (CR.GT.0.0)    GO    TO    20 

P=-2.0*CR 

C=1.C-P 

CC    100    1=1, NV 

CALL    GGE0T(CXS,1,0,WK1 ,IG) 

CALL    RANDOM( IL,U, 1) 

PN=1.0/IG 

Yd  )  =  (0=*(U**RN)  )/(1.0-(P*(U**RN)  )  ) 
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CALL    RANDCMni,G,IG  ) 

V»  =  1.0 

CO     IC    K=1,IG 

^  =  W=»G(K) 
IC    CGMINLE 

2{ I)=W**( 1,0-P) 
ICO    CCNTINUE 

n=TLRN 
20    WRIT5(6,30)CP 
30    FCRN4T(5X,»THIS    METhCO    NOT    ALLOW    CR=»»F7.?) 

RETIPN 

END 


SUBROUTINE    MAR  EXP ( CR ,NV  ,Y, Z » I  X) 


CIMENSION    Y(NV),Z(NV),E(3) 

IP  (CR.LT.0.0)     GO    TG    20 

R12=(2.C*CR)/(CR+1.C) 

Pl=1.0-R12 

R2  =  P1 

DO    5C    1=1, NV 

CALL    EXP0N(IX,E,3) 

E(1)=E(1)/R1 

E(2)  =  E(  2)/R2 

E(3  )=E(3)/R12 

>( I)=AMIN1(E(1)  ,E(3  )) 

Z(  I)=AMIN1(E(  2),F(3)) 
50    CCNTINUE 

R  =  TLRiN 
20    WRITE(6,3C)CR 
30    FCRNAT(5X,'TI-1S    METFCD    NOT    ALLOW    CP  =  »,F7.3) 

RETLRN 

END 


SUEPOUTINE    GAVEXP(CPtNV,Y,Z,IS) 

CIMCNS!0NY(NV),Z(NV),WK1(1)»WK2(1CC) 

CCLELE  PRECISION  OS E EDI ,0SEED2 , DS EEC3 

CSEED1  =  IS*0. 2  +  123456 

DSEEC2=IS*0.5n234 

CSEED3  =  IS*0. 9  +  123 

lF<CR.LT.-0.5  .OR.  CR.GT.O.OJ  GO  TC  20 

C=-2.0*CR 

P=1.C-Q 

CO  ICC  1=1, NV 

CALL  GGEOT< CSEED1,1,P,WK1, IG) 

A  =  IG 

e=p 

CALL  GGAMS(CSEED2,A,e,l  ,WK2,GA) 
2(1 )=GA 

CALL  GGUBS(CSEED3, 1,U) 
PN=1.0/IG 

XI  )=ALOG((l.C/(t>*U**RN))-(C/F)  ) 
100  CCNTINUE 
PETLRN 
20  VRITE(6,3C)CR 

30  FGPV4T{5X,»THIS  METhCD  NOT  ALLOW  CR=SF7.3) 
RETLRN 
END 
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